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PREFACE 


A  numerical  investigation  of  hurricane  induced  water  level  fluctuations 
in  Lake  Okeechobee,  Florida,  is  conducted  from  both  forecasting  and  design  per 
spectives.  The  report  presents  the  development  and  application  of  a  combined 
hurricane  and  hydrodynamic  modeling  system  complete  with  graphics.  The  hydro- 
dynamic  model  contains  both  long  and  short  wave  components. 

Project  administration  and  funding  were  provided  by  the  US  Army  Engineer 
District,  Jacksonville. 


The  numerical  investigation  was  conducted  at  the  US  Army  Engineer 
Waterways  Experiment  Station  (WES)  in  the  Coastal  Engineering  Research  Center 
(CERC)  from  June  1982  through  June  1985  under  the  general  supervision  of 
Drs.  R.  W.  Whalin  and  L.  E.  Link,  former  Chief  and  Assistant  Chief,  respec¬ 
tively,  and  Dr.  J.  R.  Houston,  former  Chief,  Research  Division  (RD)  ,  and 
now  Chief,  CERC.  This  report  was  prepared  by  Dr.  R.  A.  Schmalz,  Jr.,  RD . 

The  author  would  like  to  acknowledge  the  contributions  of  Mr.  David  L. 
Leenknecht  for  providing  the  short  wave  modeling  concepts,  as  well  as  pro¬ 
viding  preliminary  versions  of  the  graphics  system  to  be  used  in  the  Coastal 
Modeling  Package,  Mr.  Charles  C.  Curren  for  work  done  in  the  design  storm 
software,  and  Mr.  Mark  D.  Prater  for  providing  results  of  his  work  on  wind 
stress  and  drag  coefficients.  In  addition,  the  following  Jacksonville  Dis¬ 
trict  personnel  are  acknowledged  for  greatly  assisting  in  this  project. 

Mr.  Joseph  E.  Gurule,  served  as  project  manager  and  was  instrumental  in 
formulating  the  overall  study  approach  in  conjunction  with  the  technology 
transfer  process.  Mr.  Kaiser  Edmand  developed  levee  alignment,  toe  and  crest 
elevations,  and  side  slopes  as  well  as  digitizing  the  global  grid. 

Mr.  Kline  Bentley  provided  water  level  and  wind  records  for  seiches  in  August 
1969  and  1968  as  well  as  for  Hurricane  David. 

Commander  and  Director  of  WES  during  publication  of  this  report  was 


Technical  Director  was  Dr. 


Robert  W.  Whalin. 
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PART  I :  INTRODUCTION 


1.  This  report  investigates  hurricane  induced  water  level  fluctuations 
in  Lake  Okeechobee.  In  response  to  a  letter  dated  21  August  1981  from  SAJ  and 
additional  telephone  conversations,  a  meeting  was  held  on  13  October  at  WES 
with  SAJ  and  SAD  personnel  to  outline  the  scope  of  the  numerical  investiga¬ 
tion.  Based  upon  results  of  the  meeting,  WES  proposed  on  4  December  to 
develop  a  combined  hurricane  and  hydrodynamic  modeling  package  to  be  employed 
in  both  a  design  and  forecasting  mode.  Within  the  design  mode,  both  the  tra¬ 
ditional  design  hurricane  and  more  recent  joint  probability  methods  were  to  be 
incorporated.  The  hurricane  submodel  was  to  contain  both  a  direct  parametric 
and  planetary  boundary  layer  option.  The  hydrodynamic  submodel  was  to  contain 
a  short  wave  description  including  wave  runup  and  overtopping.  Upon  further 
review  by  SAD  and  SAJ,  the  project  scope  was  extended  to  consider  extreme 
hurricane  design  conditions  for  purposes  of  dam  safety  and  emergency  action 
planning . 

2.  The  above  components  of  the  study  are  presented  in  this  report 
through  the  following  structure.  Part  II  outlines  the  hurricane  model  theo¬ 
retical  development  and  numerical  implementation.  Both  direct  parametric  and 
planetary  boundary  layer  approaches  are  developed  in  a  generalized  hurricane 
model.  The  theoretical  development  of  the  long  and  short  wave  equations  are 
presented  in  Parc  III.  Detailed  numerical  implementation  of  the  hydrodynamics 
is  the  thrust  of  Part  IV.  The  seiche  calibration  and  verification  of  the  bot¬ 
tom  friction  mechanics  in  the  hydrodynamics  is  developed  in  Part  V  independent 
from  the  surface  wind  stress.  In  Part  VI,  the  wind  stress  in  the  hydrodynamics 
and  hurricane  model  are  calibrated  using  the  August  1949  hurricane.  Both  the 
planetary  boundary  layer  and  standard  project  hurricane  windfield  formulations 
are  considered.  In  Part  VII,  the  October  1950  hurricane  is  used  to  further 
substantiate  the  hydrodynamics  by  employing  the  standard  project  hurricane 
approach.  In  Part  VIII,  Hurricane  David  is  simulated  with  the  present  levee 
conf igurat ion  again  using  the  standard  project  hurricane  approach.  In  Part  IX, 
the  development  of  design  hurricanes  for  Lake  Okeechobee  is  presented.  In 
addition,  the  results  for  a  representative  breech  for  a  probable  maximum 
hurricane  are  presented  to  demonstrate  the  model's  capability  to  predict 
flooding  potential.  Joint  probability  method  design  considerations  are 
developed  in  paragraphs  275-277  in  conjunction  with  the  second  report  in  this 


series.  In  Part  X,  results  are  summarized  and  a  set  of  conclusions  drawn. 
Directions  of  future  work  by  the  District  are  indicated  in  the  context  of 
technology  transfer  throughout  the  report. 

Appendicies  A  and  B  contain  supporting  tables  and  plates,  respectivel 
Appendicies  C-G  contain  computer  code  listings  and  documentation.  All 
appendicies  appear  in  a  separate  volume. 


PART  II:  HURRICANE  MODEL  THEORY  AND  NUMERICAL  IMPLEMENTATION 


Overview 


3.  Traditionally,  the  US  Army  Corps  of  Engineers  has  employed  paramet¬ 
ric  approaches  in  determining  wind  and  pressure  fields.  A  Standard  Project 
Hurricane  (SPH)  methodology  or  direct  parametric  approach  has  evolved  from  the 
work  of  Graham  and  Nunn  [1]  and  been  most  recently  updated  by  Schwerdt  et  al. 
[2].  Recent  research  has  focused  on  numerical  modeling  of  atmospheric  circu¬ 
lation  through  parameterization  of  the  planetary  boundary  layer,  the  so-called 
PBL  approach.  The  theoretical  development  and  numerical  structure  of  these 
two  approaches  is  presented  initially  in  turn  below.  The  numerical  implemen¬ 
tation  employed  in  a  generalized  hurricane  model  is  presented  next  followed  by 
a  discussion  on  general  g  1  interpolation  and  joint  probability  method  storm 
geometry  sections. 


Direct  Parametric  Approach 

4.  Within  this  methodology  design  hurricane  conditions  have  been  devel¬ 
oped  by  considering  the  general  characteristics  of  previous  severe  hurricanes 
in  light  of  improved  knowledge  of  the  hurricane  phenomena  in  general.  Reid 
[3!  implemented  the  1972  SPH  methodology  [4]  on  a  digital  computer.  The  pro¬ 
gram  developed  allowed  the  user  to  specify  a  design  hurricane  through  input  of 
non-time  varying  hurricane  parameters  (central  pressure  depression,  radius  to 
maximum  winds,  maximum  wind  speed,  storm  track,  storm  forward  speed,  and 
azimuth  of  maximum  winds).  Alternatively,  by  time  varying  the  above  param¬ 
eters  as  well  as  specifying  an  effective  radius  to  maximum  winds,  observed 
hurricane  windfields  could  effectively  be  matched.  As  such  the  SPH  method¬ 
ology  under  the  second  option  is  descriptive  and  not  predictive  and  therefore 
is  not  suitable  in  its  current  state  for  predicting  windfields  for  hypothet¬ 
ical  hurricanes  to  be  considered  in  a  joint  probability  method  format  for 
design . 

5.  Perhaps,  the  most  widely  used  parametric  hurricane  windfield  model 
is  the  Tetra  Tech  Model  (TTM)  [5]  as  employed  for  the  coastal  flood  insurance 
•  gram  of  the  Federal  Emergency  Management  Agency.  Within  this  model,  the 
user  specifies  central  pressure  depression,  radius  to  maximum  winds,  forward 


speed  and  direction  and  landfall  point.  The  following  variables  are  con¬ 
sidered  calibration  parameters: 

a.  3  ,  Incurvature  angle  (°) 

b.  V  ,  Far  field  wind  speed  (kts) 

c.  c  ,  Reduction  factor  (maximum  gradient  to  10  m  over  water  wind 
speed)  (-) 

6.  Both  the  SPH  and  TTM  are  structured  similarly.  A  radial  profile  of 
windspeed  is  specified  along  the  radial  of  maximum  winds.  An  asymmetry  factor 
is  specified  as  an  azimuthal  relationship  to  modify  windspeed  in  the  left  and 
right  halves  of  the  storm.  A  constant  inflow  or  incurvature  angle  is  speci¬ 
fied  to  define  the  wind  direction  field. 

7.  In  order  to  develop  the  parametric  descriptions  the  following  topics 
are  presented  in  turn  below: 

a.  Cvclostrophic  and  gradient  windfield  determination 

b.  Determination  of  the  k  coefficient 

c.  Ten-meter,  10-minute  overwater  winds 
Asymmetry  Factor  (A) 

Asymmetry  Angle  (6) 

General  Windfield  Specification 

Cvclostrophic  and 

Gradient  Windfield  Determination 

8.  Consider  the  equations  of  steady  horizontal  motion  written  in  polar 
coordinates  as  follows: 


!!l+  III  _  ^9  _  I  3P 

vr  3r  V0  rae  r  9  o  a 

C.l.l) 

3V9  3va  VPVr  1  ?? 

-  +  v.  — -  +  — • —  =  -fv  -  —  — —  (b) 

r  3r  v  r39  r  r  o  r39 


where  v  ,  v.  represent  the  radial  and  tangential  velocity  components. 

9.  Consider  the  Hvdromet  pressure  distribution  developed  by  Schloemer  [6' 


as  valid  in  a  hurricane. 


P  -  P 


-R/r  Ap  (r) 


(2.1.2) 


where 


P  =  sea-level  pressure  at  distance  r  from  the  hurricane  center 
?  =  asymptotic  pressure  (29.92  in  Hg) 

R  =  storm  radius 
Pq  =  central  pressure  (in  Hg) 

This  pressure  distribution  has  concentric  isobars  with  their  centers  at 

r  =  0  .  Thus  3P/30  =  0  .  If  we  assume  that  the  velocity  distribution  has 

circular  symmetry  for  the  case  of  a  stationary  hurricane,  then  ?v  /39  , 

9vq/3  0  =  0  .  If  v^_  =  0  ,  the  second  equation  is  satisfied.  Replacing  v„ 

with  c  the  first  equation  becomes: 

C*"  1  -P 

—  +  fc  -  -  —  =  0  (2.1.3) 

r  -  3r 

«_j  i _ i  i - » 

C  C  P 
e  o 


This  equation  represents  a  balance  between  centrifugal  (C  ) ,  Coriolis  (C  ) , 

6  O 

and  pressure  gradient  force  (P)  . 

10.  If  the  curvature  effects  dominate  the  Coriolis  effects,  the 
following  equation  is  obtained  for  the  cvclostrophic  wind,  V 

2 

!c  =  i  9P  (2il 

r  -  i  r 


If  both  effects  are  important,  c  ,  represents  the  gradient  wind.  Thus  one 
obtains  the  following  relation  between  the  cvclostrophic  and  gradient  winds: 


c_ 

r 


+ 


fc 


V~ 

c 

r 


2  2 

V  -  c  =  rfc 
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(V  +  c) (V  -  c)  =  rf c 
c  c 


V 
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c 


rfc 
c  +  V 


(2.1.5) 


C  IS 


11.  The  assumption  is  made  that  the  difference  between  V  and 

c 

small  compared  with  the  quantities  themselves.  Therefore 


c 


(2.1.6) 


for  (2.1.2)  and  (2.1.4) 
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c  P  r 


Re 


-R/r 


For  V  ,  r  =  R  ;  e.g.,  the  maximum  cyclostrophic  wind  occurs  at  r  =  R  : 


V 

cx 


K(P  -  P  )1/2 
w  o 


(2.1.7) 


where  K  =  (oe) 

Returning  to  Equation  (2.1.6)  we  obtain  the  following  relationship  for  the 
maximum  gradient  windspeed,  V  : 


I  / 7  Rf 

V  =  K(P  -P  )  '  -  (2.1.8) 

gx  wo  2 

Determination  of  the  K  coefficient 

12.  The  magnitude  of  the  maximum  gradient  wind  is  dependent  not  only  on 
the  pressure  difference,  but  also  on  the  air  density  at  R  ,  which  has  been 
included  in  K  .  The  numerical  value  of  K  depends  upon  the  units  used.  For 
wind  speed  in  miles  per  hour  and  pressure  in  inches  of  mercury,  Graham  and 
Nunn  [i]  used  K  =  73  based  upon  an  air  density  at  68°F  and  a  pressure  of 
29.53  inches.  Schwerdt  et  al.  [2]  have  assumed  the  latitudinal  variation  of 


K  is  based  on  the  variation  of  the  0.99  probability  level  sea-surface  temper¬ 
ature  for  the  PMH  and  the  0.75  probability  level  for  the  SPH.  The  assumption 
is  made  that  the  air  temperature  is  equal  to  the  sea  surface  temperature.  The 
difference  in  K  at  24°N  and  45°N  is  about  4%  and  5 %  for  the  SPH  and  PMH, 
respectively. 

Ten-meter,  10-minute  overwater  winds 

13.  Observed  maximum  LO-m,  10-min  winds  (V  )  over  open  water  in  hurri¬ 
canes  of  above  average  intensity  have  been  found  to  varv  from  75-100%  of  V 

gx 

by  Myers  [7],  This  maximum  point  will  lie  on  a  circle  of  radius  R  .  In 
general , 

V  »  F(V  )  +  A  (2.3.1) 

X  gx 

where 

=  maximum  10-m,  10-min  overwater  wind  speed 
F  =  reduction  factor  to  convert  from  maximum  gradient  wind  speed  to 
10-ra,  10-min  overwater  wind  speed 
(F  =  .9  for  the  SPH 
F  *  .95  for  the  PMH) 

V  =  maximum  gradient  wind  speed  defined  in  (2.1.8) 
gx 

A  =  asymmetry  factor  resulting  from  the  forward  speed  (T)  of  the 
hurricane 

For  a  stationary  hurricane,  the  asymmetry  factor,  A  ,  is  zero.  for  a 

stationary  hurricane  we  designate  as  V  .  Knowing  ,  we  can  use  rela¬ 

tive  wind  profile  information  to  determine  10-ra,  10-min  overwater  winds  at  any 
distance  from  the  hurricane  center. 

Asymmetry  Factor  (A) 

14.  For  a  moving  hurricane,  the  asymmetry  factor  is  non-zero.  The 
general  equation  for  A  being: 


A  **  yT  cosB 


(2.3.2) 


10 


where  y  and  x  are  two  empirical  constants  dependent  upon  units,  T  is  the 
forward  speed,  and  8  is  the  angle  between  track  direction  (-)  and  the  sur¬ 
face  wing  direction  (5  ),  measured  counter-clockwise  from  9  .  Graham  and 

a 

Nunn,  [2]  used  y  ■  .5  and  x  -  1  when  the  windspeed  was  expressed  in  knots. 
Schwerdt  et  al.  [2]  suggested  y  =  1.5  and  x  -  .63  as  more  representative 
constants.  In  Table  II-l  below  the  maximum  asymmetry  added  is  contrasted. 


Table  1 1 —  1 :  Maximum  Asymmetry  In  Wind  Speed 


Forward  Speed  (kts) 
6 


Graham  and  Nunn 
3 


Schwerdt  et  al. 
4.6 


20 


10 


10.0 


50 


25 


17.6 


Asymmetry  Angle  (8) 

15.  The  asymmetry  angle  is  determined  as  shown  in  Figure  X I —  1  below. 

As  may  be  noted,  it  is  directly  dependent  upon  the  inflow  angle,  0^  ,  defined 
as  the  angle  between  the  true  wind  direction  at  a  point  and  the  tangent  to  a 
circle  about  the  storm  center  passing  through  the  point. 

16.  In  reality,  9^  ,  varies  with  radial  distance,  r  ,  outward  from 
the  storm  and  with  angle  with  respect  to  the  storm  direction.  Nunez  and  Gray 
[9]  define  the  front  quadrant  as  being  centered  around  the  track  direction  (9) 
and  note  the  boundary  layer's  inflow  angle  magnitude  decreases  from  quadrant 
to  quadrant  in  the  following  order:  right,  front,  back,  left. 

17.  In  the  SPH  model,  neither  the  variation  in  inflow  angle  with 
respect  to  orientation  about  the  storm  direction  nor  the  radial  variation  are 
considered.  In  the  first  SPH  study  by  Graham  and  Nunn  [l],  a  value  of  20° 
from  the  hurricane  center  to  the  radius  of  maximum  winds  (R) ,  a  transition 
from  20-25°  between  R  and  1.2R,  and  25°  beyond  1 . 2R  was  advocated.  Later 
studies  [8,  4]  varied  this.  Although  25°  was  used  for  1.2R  and  beyond, 
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angles  from  0  at  the  center  increasingly  linearly  to  10°  at  R  were  speci¬ 
fied.  A  transition  from  10°  to  25°  was  used  between  R  and  1.2R. 

18.  Jelesnianski  [10]  related  inflow  angle  to  maximum  surface  winds,  R 
and  pressure  drop  (P  -  P  ) .  Nomograms  can  be  constructed  at  a  given  latitude 
at  prescribed  distances  from  the  hurricane  center. 

19.  Schwerdt  et  al.  [2]  rely  on  the  results  of  Jelesnianski  [10]  with 
the  following  assumptions: 

a.  The  SPH  is  modeled  after  Jelesnianski ' s  [10]  nomogram  for  a 
Ap  =  2.08  in  Hg. 

b.  The  PMH  is  modeled  after  Jelesnianski' s  [10]  nomogram  for  a 
Ap  =  3.34  in  Hg. 

c.  Maximum  inflow  will  occur  at  a  distance  of  3R. 

d.  Inflow  will  decrease  outward  but  remain  positive  from  3R  to  the 
periphery  of  hurricane  circulation. 

e.  Inflow  will  be  constant  with  respect  to  track  direction 
orientation. 

jf.  Inflow  will  not  vary  with  storm  forward  speed. 

£.  Inflow  will  not  vary  with  latitude. 

20.  These  results  are  stated  by  Schwerdt  et  al.  [2]  to  represent  an 
improvement  over  previous  inflow  angle  criteria  from  the  following 
perspective : 

a.  Curves  of  inflow  angle  versus  R  are  now  continuous. 

b.  Maximum  inflow  is  no  longer  a  constant  but  increases  with 
increasing  R  . 

c.  Maximum  inflow  does  not  extend  to  the  periphery  of  the 
hurricane,  which  is  in  agreement  with  Chow  [11]. 

General  Windfleld  Specification 

If  we  combine  Equations  (2.9)  and  (2.10),  we  may  obtain  a  relation- 
in  a  moving  hurricane  in  the  following  form: 

V  =  F(V  )  +  1.5  T*63cos6  (2.3.3) 

x  gx 

(Note  y  »  1.5  and  x  =  .63  ,  for  expressed  in  knots)  (F  =  9.0  SPH) 

(F  »  .95  PMH) 


21. 

ship  for 
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*  V  V* 

•  .N  \ 

22.  Equation  (2.3.3)  may  be  used  to  describe  the  wind  for  r  =  R  only. 
In  other  words,  the  equation  specifies  the  windfield  at  a  radial  distance 
equal  to  the  radius  to  maximum  winds. 

23.  The  following  relationship  is  used  to  specify  the  10-m,  10-min 
overwater  winds  at  any  distance  (r)  from  the  hurricane  center: 


V  =  V  +  1.5T‘63cos  S 
s 


(2.3.4) 


where  V  is  the  windspeed  in  a  stationary  hurricane  at  radius  r  .  Relative 
winds  profiles  based  upon  observed  hurricanes  are  developed  for  0  <  r  <  R 
and  r  >  R  in  Chapter  13  in  Schwerdt  et  al.  [2].  It  is  possible  to  also 
employ  the  Hydromet  pressure  profile  in  Equation  (2.2)  to  compute  the  profile. 

24.  Tetra  Tech  based  upon  work  by  Collins  and  Viehman  [12]  assume  the 
following  radial  distribution  for  the  maximum  hurricane  windspeed. 


V  (r)  = 
max 


max  ,  R 

—  k  lo810  — 


[.5  V  f  ;  -  i) 
max  '  R  3  ' 


r  >  R 


R  <  r  >  R/3 


R/3  >  r 


(2.3.5) 


where 

V  (r)  =  maximum  windspeed  at  radius  r  (kts) 
max 

V  =  maximum  windspeed  over  all  r  ,  r  =  R  (kts) 
max 

k  =  -.15128 
c  =  3.354 
m  =  1.607 
c2  =  1.265  x  10~3 

25.  To  completely  determine  the  windfield  the  location  of  the  radial  of 
maximum  winds  with  respect  to  the  storm  direction  must  be  specified.  Schwerdt 
et  al.  [2]  allow  the  radial  of  maximum  winds  to  be  in  the  interval  (0°,  180°) 
measured  clockwise  with  respect  to  the  track  direction.  In  the  Tetra  Tech 


model,  [5]  Che  locus  of  maximum  winds  at  any  given  radial  length  r  is  a  line 
inclined  at  (90°  +  a)  to  the  direction  of  hurricane  motion.  a  is  a  user 
specified  incurvature  angle,  which  does  not  vary  with  radial  length  r  . 


Angular  Conventions 

26.  Consider  the  hurricane  windfield  setup  sketched  in  Figure  1 1  —  2 
below.  Note  from  the  construction,  we  obtain  a  +  90°  =  6+  5  =  9  +  6  ,  and 
therefore  9=6.  Since  the  cos  x  =  cos(-x),  0  may  be  measured  either 
clockwise  or  counterclockwise  from  the  radial  of  maximum  wind. 

\ 

\ 

\ 

\ 


\ 


Figure  1 1  —  2 .  Asymmetry  Angle  Relationships 

2’’.  In  the  implementation  of  the  above  procedures  the  following 
approach  outlined  in  Figure  1 1  —  3  is  employed.  Thus,  when  the  maximum  wind- 
speed  is  specified  -5.V^.  is  added  to  convert  the  maximum  windspeed  to  equiv¬ 
alent  stationary  storm  case.  This  is  necessary  since  the  asymmetry  relations 


are  to  be  applied  directly  only  to  the  stationary  hurricane.  The  azimuth  of 
maximum  winds  may  be  specified  and  is  not  equal  to  90°  +  a  .  The  9  concept 
is  used  to  specify  the  asymmetry. 

18.  In  the  current  SPH  approach,  the  inflow  angle  (a)  is  specified  as  a 

constant.  Thus  the  wind  direction,  9  ,  equals  9  +  a  +  90°  as  shown  in 

w  p 

Figure  II-4  below. 


✓  Qp 


E,  X,  AXIS 


Figure  II-4.  Determination  of  Wind  Direction 


In  component  form,  we  have  for  the  windspeed: 


V  =  V  cos(Q  +  a  +  90°)  =  (-sin©  cosn  -  sina  cos6„)  Vp  (a) 
x  p  p  p  P 


(2.3.6) 


V  =  V  sin(9  +  a  +  90°)  =  (-cos9  cosa  -  sin0„  sina)  Vp  (b) 
v  p  p  p  P 


Planetary  Boundary  Laver  Approach 


General  Theory 


29.  Chow  [II]  considered  the  following  vertically  averaged  equation  of 
motion  in  vector  form  with  respect  to  coordinates  fixed  to  the  earth: 


(3.1.1) 


with 


f 


k  x 


A 

w  =  - 


1  A  '■'d  A ,  A 

—  VP  +  V- (k„  V  w)  -  —  ; w  w 

0  H  h  1  1 


where 
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at 
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I  j 
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k 
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time  change  local  Co  the  fixed  coordinates 

two  dimensional  del  operator 
vertically  averaged  horizontal  velocity 
Coriolis  parameter 

unit  vectors  in  the  x,  y,  and  z  directions,  respectively 
Horizontal  eddy  viscosity  coefficient 
Drag  coefficient 

Depth  of  the  planetary  boundary  layer 
pressure 

mean  air  density 


It  is  assumed  that  the  vertical  advection  of  momentum,  shear  stress  at  the  top 

of  the  boundary  layer,  and  the  turning  angle  are  all  zero. 

30.  The  pressure  P  is  given  as  the  sum  of  P  ,  the  pressure  field 

c  A 

representing  a  tropical  cyclone  assumed  to  move  at  a  constant  velocity  w^  , 
and  P  a  large-scale  pressure  field  of  a  constant  gradient.  The  corre¬ 
sponding  constant  geostrophic  flow,  w  is  given  by  the  following  equation: 

f Ik  x  w  »  -  -  VP  (3.1.2) 

i  g  o 


Thus,  Equation  (3.1.1)  may  be  written  as  follows: 

dw  A  A  l  A  CD  , A , A 

—  +  f  k  x(w-w)=  -  -  VP  +  V  •  (k  V  w)  -  —  w  w  (3.1.3) 

dt  1  g  pc  H  h  1  1 

31.  Let  us  now  consider  a  moving  Cartesian  coordinate  system  (x,  y) 
whose  origin  is  located  at  the  moving  low  center.  Therefore,  we  may  obtain: 


A  A 

w  =  w  +  w 

A  +  E 

w  =  w  +  w 

g  g  C 


Substituting  relations  (3.1.4)  into  (3.1.3)  we  obtain: 


(3.1.4) 


d  A 

— (w  +w)+fik*(w-w) 
dtv  c'  1  g7 


-VP  +  V 

D  c 


[k  V(w  +  w  ) ] 

H  C 


D  I  a  ,  A 

; —  |  W  +  W  |  (w  +  w  ) 
h  C  C 


(3.1.5) 


Noting  that  =  0  ^  Vw^  =  0  ,  we  obtain  the  following  simplification: 


+  f ik  *  (w  -  wg)  = 


-  ^  VPc  +  V  •  (Ic^Vw)  “  lw  +  wc  !  (w  +  w^)  (3.1.6) 


_  _  UW  JW  ,  '  '  .  o  W 

Note:  —  =  —  +  (w  +  w  )  •  Vw  *  —  +  w  •  Vw  +  w  •  Vw 

dt  3t  c  3t  c 


sw  \ 

—  ^  +  W  •  VW 


We  observe  in  (3.1.6)  that  w  is  the  horizontal  wind  relative  to  the  low  cen¬ 
ter  and  (3v/3t)  (is  the  local  acceleration  in  the  moving  coordinate  system. 
Since  7P^  is  time  invariant  in  the  moving  coordinate  system,  the  windfield 
may  be  determined  as  a  steady  state  solution  [(3w/3t)  ♦0]  in  this  system. 

Calculations  are  actually  made  as  an  initial  value  problem  in  which  time 
marching  is  performed  until  the  windfield  becomes  approximately  steady. 

32.  Let  us  consider  w  -  (u,  v)  ,  w  =*  (u  ,  v  )  ,  and  w  =  (Q  ,  v  ) 

g  g  g  c  c  c 

then  Equation  (3.1.6)  may  be  written  in  the  following  component  form  noting: 


■f  ,  k  *  (v  -  w  )  = 
g 


-  (u  -  u  )  f 


(u  -  U  )  (V  -  V  )  0 

g  g 


(3.1.7) 


=  f  v  -  P  +  H  -  F 
dt  u  u  u 


3—  ■  fu  -  P  +  H 
dt  v  v 


l_/k  3u\  +  3_/  3u\ 

.'X  ^  H  3xy  3y  ^  H  3v  / 

-(k  4r\  +  f  (k„  Is) 

:x  I  H  ox!  dy  V  H  oy  J 


fj(k„  :?)- 

1  3F 
,  1  c 

rv  +  -  - — 

g  =  3X 
r  ,  l  _ C 

g  0  3v 


(u  +  u  )  “  +  (v  + 
c 


A  0 

V  J 


A  2  A  2 

—  (u  +  v  )  +  (v  +  v  ) 

he  c 


,  A  . 

(u  +  u  ) 
c 


,  A  . 

(v  +  V  ) 
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Note  that  the  terms  depend  on  the  wind  relative  to  the  fixed  coordi¬ 

nates.  The  non-linear  lateral  diffusion  coefficient  k u  is  expressed  fol- 

H 

lowing  Smagorinsky  [13]  as  follows: 


k  u2  /£ 
=  2k  TT 

M  o  2 


(f)2  I'D*£I 


(3.1.8) 


where 
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33.  In  practice  the  geostrophic  wind  components  are  specified  relative 

to  the  fixed  coordinates,  thus  in  the  expressions  for  P  and  P  above: 

u  v 

f  A  f  A 

U  =  U  -  U  ,V  =  V  -  V 

g  g  c  g  g  c 


where 


(  J  ,  Vg) =  the  specified  geostrophic  wind 
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34.  Caraone  [14]  has  considered  an  alternative  form  of  Equation  (3.1) 

in  which  the  turning  angle  is  non-zero.  A  revised  parameterization  of  the 
planetary  boundary  layer  as  outlined  by  Arya  [15]  is  employed  to  compute  the 
drag  coefficient  as  well  as  the  turning  angle,  3  . 

35.  The  surface  drag  relations  may  be  computed  with  respect  to  the  mear 
(vertically  integrated)  planetary  boundary  layer  wind  as  indicated  in  Equatio; 
set  (3.1.9)  below. 


•) 


C  V 
D  m 


where 


V“  =  u“ 

m  tn 


drag  coefficient 


(3.1.9) 


,  v  -  mean  wind  compounds  along  the  direction  and 
n  m  perpendicular  to  the  surface  drag 

l’“  =  'jf fiction  velocity 

A 

ku  /U.  *  -(In  z  +  A  ) 


kv  /U  =  -B  sign  f  (Note  sign  f  ,  for  our  case,  +) 


3.  k(~ 


vm 


vo 


)/5  =  -(In  z  +  c  ) 


m 


In  2  and  3  above  z  =  z  /h  where  h  is  the  height  of  the  boundarv  layer 

oo 

(assumed  spatially  uniform)  and  a  surface  roughness. 

36.  Arya  [15]  presents  the  following  relations  for  A^  ,  3^  and  : 

\ 


A  =  In  f-h/L)  +  1.5 
m 


3 

m 


1.8  exp  ( . 2h/L)  > 


h/L  ^ 


2 


fh/V*  =  1 


(3.1. 10a) 


C 


In  (-h/L)  +3.7 


(unstable ) 


- . 96 (h/L)  +  2.5 


3  =  .  J  "  h  L )  +  1 .  1 


=  - 2 . 0 ( h / L )  +  2.7 


l 3  .  1 .  1  Ob ) 


(stable) 


For  -2  <.  h/L  <  2  ,  linear  interpolation  between  unstable  conditions  at 

h/L  =  -2  and  stable  conditions  at  h/L  =2  is  empldved  to  determine  A  , 

m 

3  ,  and  C 

m  m 

37.  Tn  the  above,  h  ,  represents  the  height  of  the  boundary  layer  and 
is  considered  by  Cardone  [14]  as  a  constant.  The  Men  in-Obukhov  length,  L  , 
is  computed  as  follows: 


jvo  L>3 


(3.1.11) 


where 


(6u.) 

3  o 


=  Monin-Obukhov  length 
s  Friction  velocity 
a  von  Karman's  constant 
=  Gravity 

=  Average  temperature  at  the  surface 

=  Time  average  of  the  product  of  the  temperature  and  vertical 
velocity  component  fluctuations  at  the  surface 


2 

9  u: 

Since  9^  =  (Su^g/L:*  ,  L  =  — ^  -g -  .  For  an  unstable  atmospheric  condition, 

(9u,)  >o--?  >o->‘L<o.  If  we  consider  the  parameterization  involving 

Jo  * 

9  ,  the  '  following  relationship  for  L  is  obtained: 


A  2 

9  (In  z  +  C  )  V. 
vo _ o  m  * 

K2g  (9  -  9  ) 

vn  vo 


(3.1. ,2) 


s  -  A,  *  '  A-  »*' 


If  (In  z  +  C  )  >  o  ,  Chen  the  sign  of  L  is  determined  bv  the  sign  of 
on- 

(9  -  e  ). 

vm  vo 

38.  The  values  of  9  ,9  -  8  are  specified.  For  a  given  V 

vo  vm  vo  °  m 

and  iteration  values  for  U  and  C  ,  compute  L  and  h/L  and  z  .  The 

*  2  m  0 
roughness  length  z  =  a/g  U.  ,  where  a  is  considered  equal  to  .035.  The 

o  x 

surface  corresponds  to  open  ocean.  Then  determine  A  ,  B  .  and  C 

mm  m 

Compute  an  iterated  value  of  V  ,  V  as  follows: 

mm 


(3.1.13) 


(Note  Vu  is  determined  from  U*  ,  A  ,  B  ,  and  C  ) .  The  procedures  are 
m  ’'mm  m 

iterated  until  and  are  tolerably  close  using  the  following  inverse 

interpolation  scheme.  Let  DUV^  =  (v^  -  corresponding  to  U*i  •  If 

these  two  quantities  are  known  for  i  =  1  and  2  ,  then  an  updated  value  for 

i'  ,  u  is  calculated  with  reference  to  the  Figure  IT-5  below  as  follows. 

*  n 


Figure  II-5.  Inverse  Interpolation  Scheme 


l’n  is  limited  with  the  range  y— 5—  ,  2U^  /.  Thus  in  the  situation  shown  in 

L’s  2 

Figure  1 1- 5,  U  =  . 

The  drag  coefficient,  ,  and  turning  angle,  6  ,  are  then  determined  as 

follows : 
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o  m 


tan 


-1 


-V 
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B  >  0 


(3.1,15) 

(b) 


39.  Cardone  [14]  seeks  now  to  modify  the  third  term  (friction  term)  on 
the  right-hand  side  of  the  motion  equation  to  account  for  the  fact  that  (in 
the  northern  hemisphere)  the  surface  wind  in  a  hurricane  is  turned  inward  with 
respect  to  mean  wind  over  the  boundary  layer. 

40.  In  the  parameterization  of  the  planetary  boundary  layer  employed 

above,  a  coordinate  system  whose  x-axis  is  aligned  in  the  direction  of  the 

surface  wind  is  employed.  In  this  system,  V  <  0  as  indicated  in 
-  tn 

Figure  II-6. 
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Figure  1 1— 6 .  Planetary  Eoundary  Layer 
Parameterization  Coordinate  Svstem 


In  Equation  (3.6),  a  Cartesian  coordinate  system  whose  y-axis  is  aligned  in 
the  direction  of  North  and  whose  origin  moves  with  the  center  of  the  storm  is 
employed.  The  situation  is  shown  in  Figure  II-7  for  an  arbitrary  field  point 
in  this  system. 
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Figure  I I— 7 .  Planetary  Boundary  Layer  Computational 
Coordinate  System 

41.  In  order  to  characterize  the  conditions  at  the  surface,  V  is 

-  m 

rotated  through  6  to  obtain  Vg  .  The  friction  effect,  F  ,  is  then 


resolved  In  the  -  y^  system  as  follows: 


r  v  v 

h  s  s 


•Jui  +  v2i 


(3.1.16) 


-C  -  2 

=  —  |v  |  cos  s 
h  s 


r  D  I  2 

F  =  r—  ;  V  sin  s 
y  h  s 


where 


=  sin  (a  +  3)  =  sin  a  cos  S  +  sin  3  cos  a 
=  cos  (at  +  3)  =  cos  a  cos  3  -  sin  a  sin  6 


we  obtain: 


U1  +  V1 


sm  a  = 


cos  3  = 


2  2 
U1  +  Vi 


-(In  z  +  A  ) 

_ o  m _ 

B2  +  (In  z  +  A  ) 2 
m  o  m 


(3.1.17) 


sin  3  = 


/B2  +  (In  z.+  A)2 
'  m  o  m 


Using  the  above  formulae,  we  finally  obtain: 


cos  s  = 


sin  s  = 


-(In  z  +  A  )  U,  -  V  B 
_ o  m  1 _ l  m 

/u?  +  V2  / B2  +  (In  z  +  A  ) 2 
1  1  \j  m  o  m 
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(3.1.18) 
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B  +  (In  z'  +  A  )‘ 
»  m  o  m 


42.  In  implementing  the  turning  angle-drag  coefficient  procedures,  a 

tabular  approach  is  followed.  For  j V  |  e(.8,  80  kts)  ,  and  6  are 

calculated  for  a  water  surface.  Thus  for  an  arbitrary  field  point,  V  |  is 

m 

known  Cat  the  previous  time  step)  and  and  8  may  be  determined  from  the 

tabular  values. 


The  Posing  of  the  Initial  Value  Problem 

43.  Since  the  pressure  gradient  is  time  invariant  in  the  moving 
Cartesian  coordinate  system  fixed  to  the  storm  center,  the  local  acceleration 
in  this  system  is  zero  and  the  solution  to  the  motion  equation  is  a  steadv 
state  solution.  The  steady  state  solution  is  obtained  by  formulating  an 
initial  value  problem  and  marching  the  computations  in  time  until  the  computed 
windfield  becomes  steady.  In  this  section,  we  consider  the  specification  of 
the  initial  windfield  and  the  exterior  boundary  conditions. 

Initial  Windfield  Specification 

44.  Let  us  consider  the  equations  of  steady  horizontal  motion  in  polar 
coordinates  as  written  below: 

(a) 
(3.3.1) 

(b) 

where  v^  ,  v„  represent  the  radial  and  tangential  velocity  components. 

45.  Consider  a  circular  pressure  field  as  in  a  hurricane  (concentric 

isobars  with  their  centers  at  r  =  0)  .  Then  3 P / 3 0  =  0  .  If  we  assume  that 

the  velocity  distribution  has  circular  svmmetrv,  then  3v  / 3 a  =  3v „/99  =  0  . 

•  r  3 

If  we  take  v^_  =  0  ,  the  second  equation  is  satisfied.  Replacing  va  with  c 
the  first  equation  becomes: 


3v  3v  v  , 

_ r  _ r  9  £  1  3p 

v  -  +  v  -  -  —  »  fv  -  —  — £- 

r  3r  6  r39  r  9  3  3r 


„  Ilu  v  Ha  _  Vi 

r  3r  9  r39  r 


■fv  -  -  ^ 

r  o  r39 


—  +  fc  -  4  “  0  (3.3.2) 

r  o 

L—J  l—l  l _ l 

Ce  Co  P 

This  equation  represents  a  balance  between  centrifugal,  Coriolis,  and  pressure 
gradient  forces.  Consider  the  case  of  cyclonic  flow  around  a  low.  The  force 
balance  is  represented  as  sketched  below  in  Figure  II-8. 

|c 

t 
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Figure  I I —8  Cyclonic 
Flow  Force  3alance 


The  solutions  to  the  above  equation  are  as  follows  and  represent  the  so-callei 
gradient  wind: 


c 


(3.3.3 


For  c  =  v.  >  0  ,  a  counterclockwise  motion  (cyclonic  in  the  Northern  Hemi¬ 
sphere'!,  is  obtained.  Thus  the  plus  root  in  the  above  equation  is  considered 
to  represent  this  flow  condition.  For  a  low  3P/3r  >  0  and  the  square  root 
can  never  become  imaginary  so  that  all  values  of  pressure  gradient  are 
possible. 

46.  The  gradient  wind  condition  is  specified  as  the  initial  windfield 
condition  over  the  computational  domain.  I:  order  to  determine  the  gradient 
wind  at  an  arbitrary  field  point,  the  radial  pressure  derivative  must  be 
determined.  This  is  determined  by  assuming  a  circularly  symmetric  pressure 
field  defined  as  follows: 


P(r) 


P  +  (P  -  P  )  e  R/r 
o  no 


Therefore , 


ir 


(P  -  P  )  R 
n  o 


-R/r 
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Thus  the  gradient  wind  becomes: 
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/  2  2 

Jl  +  \1_L_ 

2  y  4 


r(p  -  ?) 

n  o 


-R/r 
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n .  3 .  - 


For  large  r  ,  this  formula  represents  c  as  the  difference  of  two  nearly 
equal  terms,  as  a  result  the  following  equivalent  form  is  used  hv  fardone 

[1M: 


29 


(3.3.3) 


Tfr  + 


R(  P  -  P  ) 

Ufr  2  +  - " - e"R/r 

2  rp 


Outermost  Boundary  Conditions 

47.  At  the  outermost  boundary,  the  acceleration  and  horizontal  diffu¬ 
sion  ire  assumed  negligible.  A  balance  between  Coriolis  force,  pressure 
I-idieot  force  and  surface  friction  is  assumed  as  shown  below. 
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V  )  - - <7P  -7—  V  +  V  (V  +  V  ) 
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(3.4.1) 


„  A  A  A 

Considering  unit  vectors  i  ,  i  ,  and  k  in  the  x,  y.  and  z  directions,  the 


above  equation  is  written  in  component  form  as  follows: 


(3. a. 2) 
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-*B.  In  order  to  develop  the  computational  method,  it  is  necessary  to 
use  the  following  relationships 
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rom  Equation  (3.4.8)  obtain  the  following  supplemental  relations: 
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(3.4.5) 
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(a) 

(3.4.7) 

(b) 

(a) 

(3.4.8) 
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(a) 


From  Equation  (3.4.8)  (a)  and  (b)  we  also  obtain: 


px  =  fv  -  r  \  u2  +  v2  u 


Py  =  -fu  -  r  J"2  +  v2  v 


(3.4.10) 


Thus,  we  may  obtain  after  some  algebra: 
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(3.4.11) 
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and  from  Equation  (3.4.11) 
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Fhus  obtain: 
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(3.4.12) 


u-  +  v-  =  (Px  +  Py)/  | -  (PX  +  V  J 


(3.4.13) 


Equation  (3.4.13)  relates  u  and  v  to  P  ,  P  ,  C  ,  h  ,  and  f  . 

x  v  D 

Let  us  consider  Equation  (3.4.8)  (a).  Multiply  by  f  and  subtract 
Cq  f~ 2  2 

-  t—  P  /u  +  v  from  both  sides  to  obtain: 
h  y  \ 
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Using  Equation  (3.4.13)  (b)  we  may  obtain: 


v  t2  +  (fs)  („2  +  v2)  -  tp,  -  ^  V777  py 


(3.4.15) 


Using  Equation  (3.4.13)  next  obtain: 
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Then  to  solve  for  u  employ  Equation  (3.4.18)  (b) ;  e.g. 
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The  numerical  algorithm  is  constructed  in  the  following  manner: 


3.  BK  =  SK 

fP  -  BK  P 

V  =  Z 

f  +  SK 
(P__  +  V  BK) 


Planetary  Boundary  Layer 
Approach  Numerical  Solution 

49.  A  numerical  solution  is  sought  to  Equation  (3.1.7).  In  order  to 
resolve  the  highly  curved  flows  a  nested  grid  system  is  employed.  Finite  dif¬ 
ference  approximations  are  employed  on  this  nested  grid  system  to  effect  the 
solution.  In  this  section  we  develop  the  nested  grid  system  conventions,  the 
general  finite  difference  approximations,  and  the  treatment  of  the  inner-mesh 
boundary  conditions. 


Nested  Grid  System  Conventions 

50.  A  nested  grid  system  of  order  five  is  employed  such  that  (A^+j  = 
2A^  ,  i  =  1  ,...4)  the  grid  spacing  doubles  as  one  proceeds  outward  from  the 
inner  mesh.  The  setup  for  the  first  two  grids  is  as  sketched  in  Figure  1 1  —  9 . 
A  similar  relationship  for  the  computational  time  step  is  employed;  e.g., 

(ac^  =  2At^  ,  i  =  1  ,...4.  The  length  of  each  side  of  grid  i  is  equal  to 

20A^  .  If  as  in  the  usual  application,  A^  =  5  km  ,  then  A^.  =  A^  •  2  = 
7[(A^  •  2)  •  2 ]  •  2 }  •  2  =  80  km  and  the  entire  grid  encompasses  (1600  km)". 
Within  the  present  computer  code,  At^  =  10A  , ;  e.g.,  for  A ^  =  5  km  then 


50  seconds. 
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51.  The  finite  difference  formulation  is  the  same  for  every  grid. 
Therefore  the  formulation  for  the  nth  grid  is  presented  with  the  subscript  n 
omitted  for  ease  of  notation.  The  approximations  for  Equation  (3.1.7)  (b)  and 
(c)  are  written  in  the  following  general  form.  (Note  a  general  variable 
i(x,  v,  t)  corresponds  to  a  finite  difference  analog  f(iA  ,  jA  ,  ~)  ,  where 
t  =  kit)  . 
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where 


P  ,  =  Components  of  the  known  pressure  gradient  force  (time  invariant) 
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Finite  difference  approximations  of  the  advective  terms 
Finite  difference  approximations  of  the  horizontal  diffusion 
of  momentum 

Surface  friction  terms 


The  treatment 

of 

the  Coriolis 

force 

is  implicit. 

Equation  set 

(4.2.  1)  above 

is  solved 
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then  Equation 
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where  FTL  =  fit 


grow 


Solving  rhe  first  equation  for  v  and  substituting  the  relation  for  v 
in  the  second  equation  one  obtains: 


- — FT-[  81  =  B2  -  FTL*  uT+1 


(4.2.3) 


Solving  for  u‘  one  obtains: 


x+1  FTL*  B2  +  B 1 

1  ~  2 
(1  +  FTL  ) 


Finally  substituting  in  the  first  relation  for  v  one  obtains 


(4.2.4) 


t+1  =  B2  -  FTL*  B1 
(1  +  FTL2) 


(4.2.4) 

b 


Equation  (4.2.4)  represents  the  forms  employed  for  computing  u  and 

v‘+l  .  In  order  to  compute  B1  and  B2 ,  (uT  ,  vT)  ,  (P^  ,  P^)  ,  (H*  ,  H*)  , 

and  (FT  ,  FT)  must  be  determined.  We  discuss  the  determination  of  these 
u  v 

terms  in  turn  below. 

Finite  Difference  Approximations 
to  the  Advective  Terms  (u  ,  V  ) 

52.  Let  us  initially  define  the  following  difference  operators: 


Cx(0i,j}  A  \  i+1/2, j  “  *i-l/2,j 


6y^i,j^  A  V  i,  j+1/2  ^i,  j-1/2 


5x(*i,j}  =y^-  ^i+1/2, j  +  1/2  ~  *i-l/2.j- 


(4.2.5) 


VT7  ri-l/2,j+l/2  ^i+l/2, j-1/ 


Then  employing  these  operators,  we  have: 
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In  the  usual  application  a  *  8  =  .5  .  The  terms  times  a  represent  upstream 
differencing  while  the  terms  times  8  correspond  to  diagonal  upstream 
differencing. 

Finite  Difference  Approximation  to 

the  Momentum  Diffusion  Terms  (H  ,  H  ) 

_ u _ v_ 

53.  Let  us  introduce  two  additional  finite  difference  operators  as 


follows : 
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he  approximations  are  then  given  as  follows: 
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et  us  expand  the  first  term  in  the  approximation  for  (H^  is  analogous 
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Expanding  further  we  obtain  the  following  relations  for  and  Dg 
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A  similar  relationship  is  obtained  for  D  with  u  and  v  exchanged. 
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Then  from  the  relations  above: 
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X  X 

X 

X  X 

I.J  l.J-1 

•  •  U  ,  V 

•  • 
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40 


-vs 


at  Che 


■V"*"?  I 


Thus  we  note  Chat  the  first  term  of  H 


following  nine  locations:  [(i,j)|i,j  t 
the  expressic 
field  points, 


the  expressions  for  H  and  H  involve 
r  u  v 


involves  u,  . 

i.J 

(-1,  0,  +  1)]. 

u.  .  and  v, 

i.J  i 


and  v  . 
i.  1 

All  other  terms 
at  these  same 


in 

nine 


Pressure  Gradient 

Force  Terms  (P  ,  P  ) 

_ u  v 

54.  The  pressure  gradient  terms  consist  of  two  components,  the  geo- 
strophic  wind  and  the  hurricane  pressure  field.  The  model  user  specifies 
through  input  the  x  and  y  components  of  the  geostrophic  wind'.  The  radial 
pressure  field  in  the  hurricane  is  determined  as  previously  developed  in  the 
section  on  the  initial  value  problem;  namely, 


3P 

_ c 

3r 


(P 

n 


-  P  )  R 
o 

2 

r 


-R/r 

e 


55.  It  is  necessary  to  compute  SP^/Sx  and  SP^/Sy  .  This  is  accom¬ 
plished  by  using  the  relationships  between  Cartesian  and  polar  coordinate 
systems  as  shown  below. 


x  =  r  cosQ 
y  =  r  sinS 


(a) 


3P 

c 


3x 


jPc 

3r  3x  3r  \r/  3r 


cos9 


3P 

c 


aP  ^  a 

ll  »  *£  / X\ 
3r  3y  3r  V  r  / 


(b) 
(4.5.1) 

(c) 


Determination  of  the 

Friction  Terms  (F  ,  F  ) 

_ u  v 

56.  The  computation  of  and  has  been  previously  presented  in 

earlier  work  on  the  turning  angle  concept.  In  the  computer  code,  the 
following  approach  is  used. 
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1.  Compute  SPH  »  ^  jVgj 

2.  Then 

SPH*  [U  (In  z  +  A  )  +  V , B  ] 
..  L  o  m  1  m 

F* - f==  = — 

\Bn,  +  (1"  ‘o  +  V 

SPH*  [V  (In  z  +  A  )  -  V  B  ] 
„  1  o  m  1  m 

y  =  n - x - 2 

J B  +  (In  z  +  A  ) 

\  m  o  m 


In  the  code  the  tangential  and  normal  drag  coefficients  are  expressed  as 
follows : 


(In  z  +  A  )  C_ 
_ o _ m  D 

J B2  +  (In  z  +  A  )2 
\  m  “ 


o  m 


CN 


B  C_ 
m  D 


I B2  +  (In  z  +  A  )2 
■*  i  m  o  m 


Ihus,  we  finally  obtain: 

1$  I 


F  = 


(U,CT  +  V,  C„) 


x  h  IT  1  N 


F  = 


( v , c_  +  U,  C  ) 


v  h  'IT  I  N 


(4.6.1) 


(a) 
(4.6.2) 

(b) 


Inner  Mesh  Boundary  Conditions 

57.  The  computational  sequence  is  constructed  such  that  one  proceeds 

from  the  outer  grid  to  the  inner  grids.  To  advance  the  solution  At,  ,  the 

o 

computations  are  performed  in  the  following  manner: 
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sequence  No. 


Grid  Nest  Number  Sequence 


The  computational  sequence  may  be  written  in  a  shorter  form  using  factorial 
notation  as  follows:  5  !  1  !  2  !  I !  3  !  1  ! 2 ! L ! U ! 1 ! 2 ! 1 ! 3 ! 1 ! 2  !  1  !  .  These  numbers  appear 
in  array  LEVEL  and  are  passed  to  Subroutine  COMQL'T  to  control  the  computa¬ 
tions.  We  note  to  advance  the  solution  At,.  :  1  6a t  ^  ,  8<it7  ,  t  ^  ,  2At7 

computations  are  performed.  The  computations  are  performed  such  that  one 
advances  in  time  from  the  innermost  to  the  outermost  grids.  In  this  manner, 
the  exterior  i  outside)  boundary  condition  for  grid  nest  i  may  be  determined 
from  the  values  on  grid  nest  i+1  .  The  interior  (.inside)  boundary  condition 
for  grid  nest  i  may  always  be  determined  or  in  fact  is  always  known,  since 
it  involves  performing  the  finite  difference  computations  extrapolating  the 
grid  corresponding  to  grid  nest  i  within  grid  nest  i-1  .  In  advancing 


computations  on  grid 


the  new  values  (L'N,  VN  arrays)  on  grid  i-1  must  be 


at  the  previous  time  level  on  grid  i  for  the  above  extrapolation  process  to 
remain  valid.  The  computational  sequence  insures  that  this  is  indeed  the 


case . 


VS.  For  grid  nest  5  the  exterior  boundary  conditions  correspond  to 
the  far  field  conditions  previously  presented. 

19.  ’.'he  exterior  boundary  conditions  for  grid  nest  i  are,  as  previ¬ 

ously  outlined,  determined  from  conditions  on  grid  nest  i+1  .  A  spatial 
interpolation  is  needed  since  the  points  on  grid  nest  i  are  twice  as  dense 
as  on  grid  i+1  .  A  temporal  interpolation  may  or  may  not  be  required.  If 
the  grid  time  level  (to  be  advanced)  on  i  is  equal  to  the  advanced  time 
level  on  grid  i+1  no  temporal  interpolation  is  required;  otherwise,  the 
boundary  values  on  grid  nest  i  are  determined  from  the  average  of  the  values 
at  the  old  and  advanced  time  levels  of  points  on  grid  nest  i+1  .  The  first 
grid  nest  number  in  each  of  the  16  factorial  sequence  presented  above  requires 
no  temporal  interpolation.  All  other  grid  nest  levels  in  each  of  the 
16  factorial  sequences  require  temporal  interpolation. 

60.  Under  the  above  approach  it  is  necessary  to  update  the  old  time 
level  arrays  (U,  V)  on  a  given  grid  nest  only  when  the  time  level  on  that  grid 
nest  is  to  be  advanced. 

61.  It  is  also  instructive  to  examine  the  range  of  indices  involved  in 
the  exterior  boundary  transfer,  interior  boundary  transfer,  and  interior 
region  computations.  The  indice  ranges  are  presented  in  Table  1 1  —  2  below. 

Table  I 1-2 

Range  of  Indices  for  Boundary  and  Interior  Regions 

I  (X  Direction)  Range  J  (Y  Direction)  Range 

Interior  Region 

Nest  1  2-20 VJ  2-20*1 

Nests  2-5  2-6  ,  16-20*J  2-6,  1 6-20 VI 

external  Boundary  1,  2I*J  1,  21*1 

Interior  Boundary 


Nest  1 
(Nests  2-5) 


N/A 

7,15  Jc (7  ,  15) 


N/A 

7,  15  It  (7 ,  15) 


Finite  Difference  Solution  Adjustment 


<•  * 


63.  If  one  considers  the  case  in  which  w  w  ,  K  ,  and  C  are  all 

'-go  D 

zero,  the  analytic  solution  corresponds  to  the  gradient  wind  field  determined 
from  the  pressure  field.  Since  there  is  no  friction  (C  =  0),  the  radial 

component  of  the  wind  is  exactly  zero. 

63.  In  the  numerical  solution,  a  negative  radial  component  of  wind  is 
present,  e.g.,  the  numerical  solution  produces  winds  with  excessive  inflow. 
Cardone  et  al.  [1A]  have  rotated  the  numerical  solution  winds  outward  by  8°  in 
order  to  remove  this  computational  bias.  This  is  accomplished  computationally 
with  reference  to  Figure  II-iO  in  the  following  manner. 


v" 


\/2  +  v2  =  V2  +  V2 

X  V  X  '  y  " 

V  .  V 

TAN  J  =  ~~  TAN  a  = 

V  X  ’  V  X 


(a) 


(b) 


J  -  >)  =  a 


Figure  11-10.  Vector  Rotation  Notation 


Consider 

sinu  =  sin  (£  -  6)  ■  sin6  cosQ  -  sin0  cosS 


(4.8.  1) 


F  rom 


N'ext 


(a)  in  Figure  TI-10 

V  =  V'  cosQ  -  V 
y  y  x 

consider , 

cosa  =  cos  (2  -  9) 


above : 
sin6 


=  cos8  cose  +  sin6 


sing 


Then 


V 

x 


I  ~  ~ 

\/V"  +  V" 
\  x  y 


V 

x 


cos9 


sinQ 


From  (a)  in  Figure  11-10  above: 
V  =  V'  cos9  +  V'  sin6 

XX  V 


(4.8.2) 


Equations  (4.8.1)  and  (4.8.2)  determine  the  coordinates  of  the  rotated  vector 

(V  ,  V  )  in  terms  of  the  numerical  solution  vector  (V*  ,  V').  Note  0=8°, 
xv  x  y 

sinr  =  .1391731  ,  cosd  =  .99026897  .  These  terms  are  found  in  Sub¬ 
routine  OL'TFLO. 


Determination  of  Surface  Stress 
and  Wind  Speed  at  Anemometer  Level 

2 

64.  The  magnitude  of  surface  stress  is  equal  to  PU^  ..where  P  is  the 

density  of  air  (dependent  on  temperature  and  pressure)  and  is  the  fric¬ 

tion  velocity  appropriate  for  the  given  surface.  The  turning  angle  associated 
with  the  mean  wind  speed  and  surface  represents  the  angle  through  which  the 
mean  wind  speed  is  turned  inward  for  the  given  surface. 

65.  The  following  logarithmic  velocity  law  is  assumed  to  hold  for  the 
given  surface  condition. 


U. 


U  *  —  In  /  — 
z  K  \  z 

o 


(4.9.1) 


where  z  =  10  meters  for  u^.the  wind  speed  at  anemometer  level. 


Numerical  Structure  of  the  Generalized  Hurricane  Submodel 

bb.  A  generalized  hurricane  model  has  been  constructed  by  incorporating 
the  parametric  approaches  within  the  five  level  nested  grid  structure  of  the 
planetary  boundary  layer  numerics.  This  approach  was  followed  in  order  to 
provide  a  direct  comparison  of  winds  calculated  at  common  grid  points.  In 
addition,  the  overall  structure  of  the  hurricane  windfields  generated  by  the 
two  approaches  may  also  be  directly  compared.  A  transfer  tape  has  been  con¬ 
structed  employing  the  following  format: 

I JPM ,  NAME,  IDAYH,  IHRH ,  NZZ,  NPOS ,  IDXN 

XLAT(I),  XLONG(I) ,  TPOS(l),  I  k  1,  NPOS) 

If  I JPM  >  0:  PINF,  CP,  RAD,  VF,  THETA,  FLINE 

For  each  snapshot:  TNZ  (N’SNAP  +  LL) ,  DELP,  RADIUS(l),  UN,  VN 

where 

IJPM  =  joint  probability  method  storm  (For  IJPM  =  0  ,  a 
standard  user  specified  storm  track  is  assumed) 

NAME  =  storm  name  (usually  storm  year  and  month) 

IDAYH  =  storm  day 

IHRH  =  storm  hour  initial  time  in  GMT 
NZZ  -  number  of  snapshots 

NPOS  -  number  of  storm  positions 

IDXN  -  minimum  spacing  of  the  innermost  grid  in  km 

XLAT(I)  -  latitude  of  storm  position  I  in  degrees  N 

XLONG(I)  -  longitude  of  storm  position  I  in  degrees  W 

TPOS(I)  -  time  in  hour  relative  to  initial  storm  hour  (IHRH) 
of  position  I 

PINF  =  far  field  storm  pressure  (mb) 

CP  =  central  pressure  (mb) 

RAD  =  radius  (nm) 

VF  =  forward  speed  (kts) 

THETA  =  line  angle  (degrees  clockwise  from  North) 

FLINE  -  storm  path  number 

TNZ  (NSNAP  +  LL)  -  snapshot  time  in  hours  relative  to  IHRH 


DELP  =  Central  Pressure  Depression  (in  Hg) 

UN  =  x-velocity  field  components  over  the  complete  five 
level  nested  grid 

VN  =  y-velocity  field  components  over  the  complete  five 
level  nested  grid 

67.  The  transfer  tape  is  accessed  by  the  Program  HGRAPH,  the  hurricane 
model  graphics  package  as  well  as  Program  LAKE,  the  hydrodynamic  model.  The 
generalized  hurricane  model  employs  a  snapshot  philosophy.  In  hindcasting  or 
forecasting  mode,  the  model  user  specifies  storm  positions  over  time  and  at 
appropriate  points  in  time  the  following  hurricane  characteristics: 

Forward  Speed  (kts) 

Direction  (Clockwise  from  North  in  degrees) 

Radius  (nm) 

Far  field  pressure  (mb) 

Central  pressure  (mb) 

For  PBL  (Steering  flow  (m/sec) 

Formulation  Steering  flow  direction  (Counter-clockwise  from  East  in 
only  degrees) 

68.  The  model  then  computes  UN  and  VN  over  the  nested  grid  for  each 

snapshot.  In  this  approach,  if  storm  characteristics  are  invariant  with 

respect  to  time  only  a  single  snapshot  need  be  computed.  The  model  user  may 
select  through  variable  IWIND  either  a  planetary  boundary  layer  windfield 
determination  or  a  parametric  approach.  Within  the  parametric  approach  either 
the  SPH  or  TTM  windfield  formulation  may  be  selected  by  specifying  the 
appropriate  value  of  IFORM.  Sample  input  data  sets  are  presented  in 

Tables  1 1  —  3  and  II-4,  respectively  for  each  approach.  If  the  model  user 
specifies  IJPM  >  0  ,  a  straight  line  storm  track  is  determined  by  the  model 
as  discussed  subsequently  in  the  joint  probability  method  storm  geometry  sec¬ 
tion.  A  test  data  set  for  joint  probability  method  storm  geometry  is  given  in 
Table  1 1  —  5 .  Within  the  joint  probability  mode  either  windfield  approach  may 
be  selected. 


Table  IT-3 

SPH  Parametric  SPH  Formulation  Input  Data  Set 


I-snE ! 

wide i ;p«-e , 

«•  *»D 

: 3<9  e«8 

36. 

3 


as .  i 

’7.57 

8 

as. 3 

”  •  8 

36.8 

81.1 

:  3 

37.85 

88.617 

'.6 

37.17 

88.817 

«  7 

37 . 37S 

81 .8583 

18 

38  .  3 

33.3 

35 

39.3 

S 

83.55 

31 

8. 

13  . 

17. 

8. 

33.5 

8. 

33  .  5 

8. 

33.5 

8. 

33.5 

8  . 

33.5 

8  •  8 
SNARE  4 

A2NX-3S. .ALP-18.  , 
8  END 
INANE  4 

AZNX-48. .ALP-18. . 
(END 
INANE  4 

AZNX-48. .ALP-11. , 
•  END 

INANE! 


3S. 

31. 

14.1 

1815. 

975. 

395 

14  . 

18 1 S . 

95  4  . 

388 

13.1 

1815. 

961  . 

388 

11.9 

1815. 

971  . 

389 

11.4 

1815. 

984. 
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AZNX-49. .ALP-18., 
•  END 

INANE 4 

A2NX-83. .ALP-18.. 
■END 


Table  II-4 

Planetary  Boundary  Layer  PBL  Formulation  Input  Data  Set 


•  NAME  1 

I B  =  1 ,  IWIND*0,  I J  PM  *  0 , 
I  END 

1349  003 


Co.  11. 

3 


35.  i 

*7  r  ** 

«  -/  / 

0. 

77.3 

1  . 

C  o  .  3 

SO.  1 

13. 

;~.o: 

30.617 

16. 

3~.  l ~ 

30.817 

17. 

:  a .  3 "  z 

31.0333 

18. 

.3 . 3 

83.3 

C" 

W  . 

Z  9 . 3 

3  2  .  -j  o 

31  . 

y  m 

1  j  . 

1  "  . 

"4  f 
w  ->  • 

3  1  . 

“1 

4  -J  -J  a 

-x  —  C- 

—  w  • 

14.1 

1013. 

J~3  . 

~  9  3 

** 

1  \  m  . 

1  4  . 

10  13. 

934  . 

3  3  9 

7  . 

143. 

-  -  f 

H  W  •  J 

13.1 

1013. 

961  . 

308 

141  . 

-1  ~t  r 

MW  •  J 

11.9 

1013. 

971  . 

33  9 

*3 

107. 

-*  -*  cr 

10.4 

1013. 

984  . 
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Table  1 1  —  5 


Toint  Probability  Method  Storm  Geometrv  Test  Data  Set 


♦NArtEI 

IB* 1 ,  IW  INO*  I  .  I  JPM=  1  , 

♦  END 

26.92  80.83 

31.88  25.2  8  1.  2  ;  . 

15.  45.  0 .  90. 

1  0  0 

♦NAnE4 

-2MX  =  5  5 . , ALP  =10.  , 

♦  End 


29.  82.5  26.5 

28.2  30.625 

29.91  28.73  12. 


69.  Due  to  the  possibility  of  temporal  variations  of  the  land/water 
boundary  of  Lake  Okeechobee  over  the  course  of  a  hurricane  event,  it  is  neces¬ 
sary  to  employ  a  dynamic  reduction  factor  methodology.  The  changing  land 
surface  will  introduce  a  time  varying  reduction  to  overwater  winds  in  the 
vicinity  of  the  boundary.  To  incorporate  this  effect  dynamically,  the  wind 
reduction  scheme  employed  by  Reid  [3]  had  been  incorporated  directly  in  the 
hydrodynamic  model.  This  approach  is  in  contrast  to  the  previous  approach 
emploved  by  Reid  [3]  in  determining  the  reduction  factors  within  the  hurricane 
model . 


Generalized  Hurricane  Submodel  Interpolation  (Grid)  Concepts 

70.  Within  either  the  parametric  or  planetary  boundary  layer  approach 
of  the  generalized  hurricane  model,  computations  are  performed  on  a  nested 
five  level  grid  in  order  to  provide  for  direct  windfield  comparisons.  Several 
levels  of  interpolation  are  performed  in  the  hydrodynamic  model  as  indicated 
in  Table  1 1 —  6.  Each  interpolation  is  presented  in  turn  below. 


•S' —  '^f\  ,*.A  .Wt.'.  .X  H  .  *7-  *'•  .9*-  .  -  ,.  ..  ft  w'W  ^  M  21, 
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Table  II-6 

nterpolation  Structure  Employed  in  the  Hydrodynamic  Model 


Temporal  Snapshot  Interpolation 
Uniform  Hurricane  Grid  Spatial  Interpolation 
Stretched  Hydrodynamic  Grid  Spatial  Interpolation 


Temporal  Snapshot  Interpolation 


71.  The  storm  position  in  latitude  and  longitude  coordinates  over  time 
with  respect  to  the  start  of  the  storm  is  written  on  the  transfer  tape. 

Within  the  hydrodynamic  model  all  times  are  also  with  respect  to  the  start  of 
the  storm.  A  temporal  interpolation  is  performed  every  KSI  time  steps  within 
the  hydrodynamic  model  to  determine  the  snapshot  windfield  over  the  five  level 
nested  grid  system  employed  in  the  hurricane  model. 


'V 


ilVA  *•  **  A  JV  -  -  k  -  J  *  -  -  k  «  -  "  «  “  j  ” 


-•  *  •  "  -  \ 


Tf'.r 


Uniform  Hurricane 

Grid  Spatial  Interpolation 


7?.  In  order  to  interpolate  the  windfield  on  the  five  level  nested  grid 
to  a  uniform  hurricane  grid,  it  is  first  necessary  to  determine  the  new  posi¬ 
tion  in  latitude  and  longitude  coordinates  of  the  storm  center.  Latitude  and 
longitude  coordinates  at  the  center  of  each  hurricane  cell  are  interpolated  at 
every  time  step  from  the  transferred  position  coordinate  history. 

73.  Consider  the  following  relations  available  from  spherical  trigo¬ 
nometry  for  the  spherical  triangle  shown  in  Figure  1 1 —  1 L  below: 


Figure  1 1—  II.  Spherical  Triangle 


Vfa+c-b)  .  (a  +  b  -  c) 
sm  2 -  sin  - - - 

(a  +  b  +  c)  ,  (b  +  c  -  a) 

sin  - = -  sin  - t - 


(6.2.1) 


hav  c  =  hav  (a  -  b)  +  sin  a  sin  b  hav  c 


(6.2.2) 


74.  We  are  concerned  with  two  points,  PO  ,  the  eye  of  the  hurricane, 
and,  PI  ,  an  arbitrary  field  point  corresponding  to  the  center  of  each  uni¬ 
form  grid  cell  at  which  wind  speed  and  direction  are  desired.  Both  points  are 
specified  in  latitude,  longitude  conventions  as  follows: 


Ut  ,  ,  PI  =  [X2  ,  LJ 


75.  As  such  C  corresponds  to  the  North  Pole  for  PO  and  PI  in  the 
northern  hemisphere  (our  region  of  concern) ,  B  corresponds  to  PI  ,  and  A 
to  PO  .  Furthermore,  we  observe: 


b  =  c°L(*1)  =90  -  A] 


a  =  cor  (a,)  =  90°  -  X, 


7  =  DLq(L1  ,  L2)  -  L2  -  (c) 


Thus  for  che  case  of  concern.  Equations  (6.2.1)  and  6.2.2)  above  become: 


,  a  sin(90°  -  a2  +  C  -  90°  +  ^/i/2  sin(180°  -  A[  -  X2  -  0/2 
tan  2  =  sin(  180°  -  \l  -  \.y  +  C) /  2  sin (90°  -  ^  +  C  -  90°  +  *,)/2 


(6.2.3) 


cos(A^  +  X2  +  C)/2 


(6.2.1) 


sin(A^  -  X2  +  C )  j 2  sin(90°  -  A^  +  X2  +  C)  jz 
sin (90°  -  (Xi  +  X2  -  C)/2  sin(X2  -  ^  +  cTj 2 


cos(A^  +  X2  -  C)/2 


(Note:  0°  <  a  <  180°  and  thus  0°  <  —  <  90°) 


hav  c  =  hav(X^  +  \?)  +  cosA^  cosX?  hav  (L2  -  L^)  (6.2.2) 


«'ote:  hav(x  -  v)  =  hav(y  -  x);  since  hav  x  =  sin  -  =  (1  -  cosx)/2  and  if 


=  hav(x)  ,  x£ (0, 180°) 


~b.  In  computation,  one  uses  Equation  (6.2.2)  to  determine  c  .  Equa¬ 
tion  •  o . _ . L  >  is  then  used  to  determine  1  . 

I".  If  c  as  computed  from  Equation  (6.2.2)  is  very  small,  then 
■  ,  =  %  and  Equation  (6.2.1)  is  not  used.  Instead,  one  considers  plane 
trigonometric  relationships. 

*8.  Consider  the  following  diagram  shown  in  Figure  11-12  in  which  the 
bearing  from  P0  to  PI  is  measured  clockwise  from  North;  the  locations  of 
?I  in  each  of  four  quadrants  is  shown  as  well  as  the  signs  of  the  coordinate 
distances.  The  horizontal  distance  (signed),  H  ,  and  vertical  distance 
'signed),  V  ,  are  given  by  the  following  relations: 


/'s''  j 


i 


PI  (X,.L,) 


‘1-  L,» 


Figure  11-12.  Plane  Triangle 

H  =  (L ^  -  L2)  cos ( X  ^  +  X0) / 2  =  (L^  -  L^)  cosXj  =  (L^  -  L2)  cosX9  (a) 

(6.2.4) 

V  =  (X,  -  Xj)  (b) 

Then 

^  =  AT AN 2  (H,  V)  . 

79.  Once  the  distance  and  bearing  are  known  for  uniform  grid  cell  cen¬ 
ter,  an  additional  procedure  is  employed  to  locate  the  four  surrounding  points 
on  the  grid  nest.  A  bilinear  interpolation  as  given  by  the  following  rela¬ 
tionship  is  used  to  interpolate  winds  to  the  uniform  grid. 


YW  ■  f  UNI  +  f  UNI 

i,j  11  II+l,JJ+l,IZ  r12  L  II+1,JJ+2,IZ 


where 


+  f21  UNIii+2,jj+i>iz  +  f 22  UNIII+2,JJ+2,IZ 


XW  =  -f  VNI  -  f  VNI 

i , j  11  II+l.JJ+l.IZ  12  II+1,JJ+2,IZ 


f  VNI  -  f  VNI 

21  II+2.JJ+1.IZ  22  II+2 , JJ+2 , IZ 


YW  .  =  y  component  of  wind  on  the  uniform  grid  at  cell 
center  i,j 

XW^  .  5  x  component  of  wind  on  the  uniform  grid  at  cell 
center  i,j 


(6.2.5) 
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f  p  >  f  O  1  »  f  ' 


=  interpolation  parameters  which  sum  to  one  and  are 
dependent  on  the  relation  of  the  cell  center  (i,j)  to 
each  of  the  four  surrounding  cells  on  the  nested 
grid;  e.g.,  (II+l,  JJ+1 ,  IZ) ,  (II+l,  JJ+2,  IZ) , 

( I 1+2 ,  JJ+1,  IZ),  and  (II+2,  JJ+2,  IZ) ,  where  II,  JJ 
correspond  to  x  and  y  cell  indices  and  IZ  is  the 
nest  number 


UNI 


II.JJ.IZ 


=  x-wind  component  at  x  index  II,  y  index  JJ,  and  nest 
index  IZ 


VNI 


II, JJ.IZ 


y-wind  component  at  x  index  II,  y  index  JJ,  and  nest 
index  IZ 


The  uniform  grid  ,is  in  the  WIFM  convention  (y  to  the  right,  x  pointed  down) 
and  as  a  result,  the  y  WIFM  axis  corresponds  to  the  x  Cartesian  coordi¬ 
nate,  while  the  x  WIFM  axis  corresponds  to  minus  the  y  Cartesian  coordi¬ 
nate,  This  coordinate  transformation  is  embodied  in  the  above  two  interpola¬ 
tion  equations. 

80.  The  pressure  -anomaly  in  feet  of  water  is  computed  using  Equation 
(2.1.2)  at  the  center  of  each  cell  of  the  uniform  grid.  This  uniform  grid  is 
employed  as  an  intermediate  grid  in  order  to  determine  a  wind  reduction  factor 
for  each  octant  of  wind  direction  for  each  grid  cell.  In  addition,  a  fetch 
associated  with  the  appropriate  octant  of  wind  direction  is  also  determined 
for  the  wave  calculations.  The  details  of  these  procedures  follow  the 
approach  by  Reid  [3]  based  upon  results  in  Reference  4.  Wind  reduction 
factors  and  fetches  are  computed  every  KSI  time  steps. 


Sketched  Hydrodynamic 
Grid  Spatial  Interpolation 

81.  The  bilinear  interpolation  procedure  previously  employed  is  used  to 
interpolate  the  wind  and  pressure  field  to  the  nonuniform  stretched  grid. 
Pressure  and  wind  velocity  components  are  determined  at  cell  centers  of  the 
stretched  grid.  The  wind  velocity  components  are  transformed  to  the  corre¬ 
sponding  stress  components  via  the  drag  coefficient  formulation  presented  in 
Part  IV. 


Joint  Probability  Method  Storm  Track  Geometry 

82.  Within  the  joint  probability  method,  straight  line  storm  tracks  are 
assumed.  Although  there  is  some  evidence,  that  parabolic  tracks  might  be 
employed  in  order  to  fit  recurving  storms,  this  approach  was  not  followed  in 
the  present  study.  Rather  the  geometry  as  shown  in  Figure  11-13  was  assumed. 


Figure  11-13.  Joint  Probability  Method  Storm  Track  e^f- 

83.  Note  the  storm  direction,  6  »  as  measured  clockwise  ••  - 

equal  to  270°  +  g  .  The  storm  direction,  a  ,  as  measured  or  r  >•  • 

'x 

from  East  (the  x-axis)  is  equal  to  180°  -  a  . 
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ru  '  simulation  length  (hrs) 
R  :  storm  radius  (nm) 

V_  ;  storm  speed  ( k  t  s ) 


Ihe  storm  position  is  then  determined  by  using  the  following  relation: 


V.  =  Y  +  (nR  cos?  -  3R  sin?)/60  +  iV_  sin(180  -  ?)/60 
10  r 


X.  =  X  +  (nR  sin?  +  3R  cos0)/(bO  cosY,) 

10  i 

-  iV_  cos(180  -  ?)/(60  cos  Y.) 
F  i 


i  =  1 , . . .nT  ,  n  £ (-2 ,  -1 ,  0,  I ,  2) 


where 


=  longitude  of  the  storm  position  at  storm  hour  i 
£  longitude  of  the  lake  center 

H  latitude  of  the  storm  position  at  storm  hour  i 
£  latitude  of  the  lake  center 
H  radius  of  the  storm  (nm) 

=  bearing  of  the  storm  line  (clockwise  from  North  in  degrees) 
=  storm  speed  (in  kts) 

=  integer  signifying  storm  hour 
3  path  number  of  the  given  storm 


These  relationships  were  subsequently  slightly  modified  as  discussed  in 
Part  IX. 

In  order  to  consider  the  influence  of  overland  storm  filling  the 
coastline  configuration  is  first  specified  using  the  following  ordered  pairs 


L?i>  t 9  ,  i  r  9  _ ,  r]  j  *  1  >  n 


-  bearing  clockwise  from  North  (°)  for  point  j 

J 

r  :  distance  from  the  shoreline  to  the  center  of  the  lake  (nm)  for 
point  j 

In  application  to  Lake  Okeechobee,  the  points  shown  in  Table  I 1—7  are  used. 

89.  A  dummy  point  (0°,  9999  nm)  or  (360°,  9999  nm)  is  added  to  the  array 
of  coastal  points  in  order  to  allow  for  determining  time  over  land  for  some 
storm  tracks.  The  actual  coastline  is  specified  between  the  angles  8.08  and 
306.47  degrees  clockwise  from  North. 

90.  For  each  storm  hour,  i  the  distance  D  and  bearing  clockwise 

from  North,  3^  ,  of  the  storm  position  with  respect  to  the  lake  center  is 

determined.  For  3 .  <  3 .  <  9  ,  a  test  is  made  on  r  =■  r .  +  (r..,  -  r.) 

3  i  J+1  J  J  +  l  J 

( i ,  -  9.)/  (9.  -  9.)  ,  the  interpolated  distance  to  the  shoreline  for  8 

l  j  j+L  j  r  l 

If  £  r  ,  then  the  storm  is  considered  to  be  overland  for  the  given  hour. 

For  each  storm  path,  the  number  of  hours  the  storm  is  overland  is  determined. 
The  total  number  of  hours  the  storm  is  overland  is  then  used  to  reduce  the 
storm  central  pressure.  An  additional  snapshot  corresponding  to  the  reduced 
central  pressure  is  then  computed  employing  the  following  relation. 

AP'  =  Ap  ( 1  -  ,013h) “ 

where 

Ap'  =  reduced  central  pressure  (in  Hg) 

IP  -  original  overwater  central  pressure  (in  Hg) 
h  =  number  of  hours  storm  is  over  land 


Florida  Coastline  Specification  With  Respect  To  Lake  Okeechobee 


Angle 

(CLW  N) 

360.00 

Segment  No.  1 

Distance 

(MM) 

9999.00 

Angle 

(clw  :;) 

306.47 

Segment  No.  2 

Distance 

(NM) 

110.02 

Angle 

(CLW  N) 

246.11 

Segment  No.  3 

Distance 

(NM) 

61.66 

Angle 

(CLW  N) 

185.11 

Segment  No.  4 

Distance 

(NM) 

103.61 

Angle 

(CLW  N) 

122.62 

Segment  No.  5 

Distance 

(NM) 

46.54 

Angle 

(CLW  N) 

8.03 

Segment  No.  6 

Distance 

(NM) 

77.57 

Segment  No.  7 

Angle  (CLW  N) 


0.00 


Distance  (NM) 


9999.00 


PART  III:  HYDRODYNAMIC  MODEL:  THEORETICAL  DEVELOPMENT 


91.  The  hydrodynamic  equations  are  developed  in  turn  for  both  the  long 
and  short  wave  model  components.  The  long  wave  equations  are  presented  in  two 
dimensional  vertically  integrated  form.  The  time  averaging  concepts  necessary 
in  describing  the  turbulence  are  not  considered  in  the  derivations. 

Long  wave  equations 

92.  The  general  equations  of  the  classical  hydrodynamics  for  incom¬ 
pressible  flow  are  given  following  Lai's  development  as  follows  [17]. 


3u  3v  3v 
lx  3y  3z 


0 


(1.1) 


0 


Du 

Dt 


PF  -  +  uA2u 

x  3x 


(1.2) 


0  *  OF  -  +  u£2v  (1.3) 

Dt  y  3y 

o  -  P F  -  +  uA2w  (1.3) 

Dt  z  3z 


with 


D_ 

Dt 


v 


3  3 

37  +  W  37 


+ 


3 


2 


*7- 


where 


x ,  y ,  z 

u,v,w 


p 

u 

p 

t 


Cartesian  coordinates 

Velocity  components  in  the  x,  y,  and  z  directions, 
respectively 

Body  forces  in  the  x,  y,  and  z  directions,  respectively 

Fluid  density 

Fluid  viscosity 

Pressure 

Time 
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EnpLoy  the  kinematic  boundary  condition;  namely,  for  F(x,y,z,t)  =0  as  a 
boundary  surface  assume  any  particle  on  the  surface  remains  on  it  implies 


21  21 
DC  =  3 1 


12 

U  _  +  V  _  +  w  . 

.’X  i  y  3  z 


0 


94.  Consider  z  =  -(x,y,t)  at  the  free  surface,  then 
F  =  z  -  "(x,v,t)  =  0  .  Hence 


DF 

Dt 


a  ,  3n  .  3 " 

0  =>  - —  +  u  —  +  v  —  =  w 
3  t  3  X  3  V 


(1.1.4) 


At  the  bottom  z  =  Z  (x,y)  . 

b 

Hence 


DF 

Dt 


3Z. 

0  =>  u  — ^  +  v 
3x 


w 


(1.1.5) 


V. 


Returning  to  (1,1.1) 

n 

/  u 


3Zt 

dz  +  u(V  ST 


u(n)  t —  +  — 
3x  jy 


f 


v  dz  +  V(Zb)  li 


3n 


v(n)  +  w(n)  -  wCZ^)  = 


(1.1.6) 


The  sum  of  the  -  terms  is  zero  from  the  bottom  boundary  condition.  The 

sum  of  the  _  terms  is  equal  to  3n/3t  from  the  free  surface  boundary 

conditions.  Thus  we  obtain: 


Denoting 


3_n  3_ 

3 1  3x 


; 

zb 


u  dz  + 


3_ 

3y 


;  I 

I 


v  dz 


0 


u  dz  and  v 


"b 


(1.1.7) 


one  obtains 
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0 


(1.1.8) 


+  -i—  [("  -  z.)u]  +  r-  [(".  -  Z  )v]  = 
?t  b  jv  b 


letting  h  =  "  -  Z,  and  dropping  the  bar  notation  with  the  understanding 
b 

henceforth  we  are  considering  vertically  averaged  quantities  one  obtains 


w+k  <hu)  (hv>  ' 0 


(1.1.9) 


Equations  of  motion 

2 

95.  Consider  the  z  motion  equation  in  which  Dw/Dt  =  0  ,  ul“w  =  0 
from  assumptions  2,  3  and  6,  respectively.  Thus  we  obtain 

OF  -  -  0  (1.2.1) 

2  3  z 


is  replaced  by  -g  due  to  the  following  assumptions 

a.  The  vertical  component  of  the  Coriolis  force  is  negligible  with 
respect  to  g  . 

b.  The  vertical  tide  generating  force  component  is  negligible  with 
respect  to  g  . 

If  we  integrate  the  above  relation  from  an  arbitrary  level  z  to  the  water 
surface  obtain 


p  (z) 


n 


p(n) 


og  dr 


(1.2.2) 


To  depth  we  consider  the  following  relations  to  hold,  in  which  the  bar  quanti¬ 
ties  are  depth  averages  and  the  prime  quantities  the  local  fluctuation  from 
these  averages. 


o(z)  =  o  +  o'  (z)  (1.2.3) 

p(z)=p+p'(z)  (1.2.4) 


(  o'  dz  = 


(  p'  dz  =  0 


(1.2.5) 


Thus  we  finally  obtain  for  (T)  the  following  expression: 


t-  (p  +  p’)  dz  =  (r  -  z  )  +  p'  t - p'  ~ 

3  X  h  3  X  r  3X  r  ;x 

Zb 


If  we  let  h  =  r,  -  z  and  assume  p'  3z  /3x  ,  p'  3n/3x  =  0  ,  then 


obtain 


Zb  b 


T7  (P  +  P')  dz  =  h 


Evaluating  (2)  note  3p(n)/3x  is  not  a  function  of  depth;  thus  we  obtain 


3o(n)  _  oo  (n)  ,  s  L  pa 

— -  dz  =  — - —  (n  -  z,  )  =  h  "r — 

3x  3x  b  3x 


where 


Pa  =  p(-.) 


Next  consider  the  iterated  integral  expression  for  (3)  as  follows: 


77  /  cg  dr  =  §  h  [p(r>  - Z)I 


3o  N  .  t  9(n  -  z) 

g  _  (n  .  z)  +  p  _ - 


(  / 

J  8  77  ' 


z)  dz  "  8  37 


,  10  h”  , , 

S  ds  =  8  37  2~  ( 


s 


2 


ds  =  -dz 

Observe  from  Leibniz  rule 


/  ^ 


3  (  '  -  z ) 


dz  =  go 


4  / 


("  -  z)  dz 


+  (r'  -  V  TT  -  (n  A  £ 


(1.2.15) 


Letting  h  =  n  -  obtain: 


'(r  -  z) 


-  3  ('•'V,  5Z» 

80  3j(-rh~ 


Ah  ^  SZb\  -  3n 

80h  fe+  177"  80h 


(1.2.16) 


and  the  evaluation  (3)  is  complete.  Assembling  all  our  results  we  obtain 


f inallv 


„  —  3p  -  ,2 

,  ?p  a  ,  -  ,  3n  .  h  3o 

h  at  =  IT  +  pgh  3^  +  8  T  3? 


Analogously,  the  expression  for  the  y  gradient  is  given  by 


3p  ,  ,  —  ,  3n  ,  h^  3o 


(1.2.17) 


,  3p  u  a  ,  ~  ,  3 n  n  30 
h  T^-sh- —  +  ogh  -r —  +  g  —  7— 

3y  3x  6  3y  6  2  3y 


(1.2.18) 


Thus  we  have  employed  the  z  motion  equation  to  evaluate  the  horizontal 
pressure  gradients  in  the  x  and  y  motion  equations.  The  expressions 
obtained  above  for  these  gradients  include  the  atmospheric  pressure  anomaly, 
the  water  surface  elevation  gradient,  and  the  density  gradient.  In  the 
present  studv,  density  gradients  will  not  be  considered. 


96.  Let  us  next  consider  the  material  derivative  (left  hand  side)  of 
0.2)  the  x  motion  equation.  Expanding  the  material  derivative  and  adding 
(1.1)  we  obta in 


.-u  ,  ?u  ,  r’U  ,  ?u  ,  /  dU  ,  3v 

•7—  +  u  -t—  ■*-  v  —  +  w  —  +  u  —  +  7— 

3t  dX  ,’v  3z  \  ox  3v 


3u  +  3  (u~)  +  3  (uv)  +  9(uw) 
3t  3x  3v  3z 


3(ou)  +  3  (  pu“)  3 ( puv)  +  3(puw) 
3t  3x  3y  3z 


(1.2.  19) 


Integrating  the  last  result  over  the  vertical  (which  also  holds  for 
compressible  flow) 


3  (ou) 


3(pu2) 


3 (puv) 


3 (puw) 


(1.2.20) 


Again  employing  Leibniz  rule 


z  n 

—  (pu)  dz  =  j  —  (ou)  dz 


+  pu  ~  -  (ou)  JQ- 
n  z, 

D 


(1.2.21) 


I  :uu ) 


X  cuu) 
& 


.  ^  5n  b 

+  CUU  -T-  -  ouu  -rr— 

2X  3x 


(1.2.22) 


( cuv)  dz  = 


3(  Cuv) 


3n  b 

+  cuv  —  -  OUV  — - 

I  By  3v 

h  Z. 

I  b 


(1.2.23) 


( Ouw)  dz 


(1.2.24) 


97.  Notice  in  the  above  3Z,  /  3t  =  0  ,  since  the  bottom  is  assumed 

b 

rigid.  Denoting  terms  in  (1.1.4)  by  i£>  ,  and  those  in  (1.1.3)  by  ^  ,  we 
obtain  the  following 


Tt  j  (“u)  dz  '  0u  J£  +  h  f  (ouu)  dz 

z  v'  z 

b  % 


©  © 

1)7  S 

.  „  b  3n  .  3 

+  Cuu  ■  - —  -  ouu  —  +  — 

7  ?x  3x  3y  > 


(duv)  dz 


.  n  b  3n  , 

+  DUV  T -  -  OUV  —  4-  OUW  -  OUW 

Zb  Y  n  3y  n  2b 


(1.2.25) 
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■W  » I1  *  M  ■  1  ■  "  ■  ' 


K-  •' 


Fv 


~ —  (uhl  c  — 


-•IX 


u  u  [1  +  2u'(z)  +  u'(z)  ]  dz 


+  “  9v 


u  v  [I  +  u' (z)  +  v'(z)  +  u'(z)v'(z)]  dz 


(1.2.31) 


Let 


j  =  _ 
h 


u  +  [u'Cz)]*-}  dz  =  i- 


J 

z. 


n 

r 


[1  +  u' (z)v' (z) ]  dz 


Noting  the  depth  integral  of  the  product  of  a  bar  and  primed  quantity  is  zero, 
we  finally  obtain 


f-  (uh)  +  ( 6uuh)  +  I—  ( Buvh) 

dt  dX  ay 


(1.2.32) 


Analogously,  one  obtains  the  following  expression  for  the  left  hand  side  of 
the  y-motion  equation  (1.3)  after  integration  over  depth 


—  (vh)  +  -I-  (8uvh)  +  (wh) 
at  3x  dy 


(1.2.33) 


where 


l 

h 


{1  +  [v'(z)P}  =  i- 


n 

1  f 

h  J 
Z, 


[1  +  u' (z)v' (z) ]  dz 


Thus  3  in  (1.2.32)  and  in  (1.2.33)  is  the  same  quantity  and  is  usually 
assumed  equal  to  unity. 

99.  Let  us  now  consider  the  right  hand  side  of  the  x  and  v  motion 

equations.  Considering  the  F  and  F  terms,  we  obtain 

x  y 


F  =  f.v  +  g  +  G  F  =*  -f.u  +  g  +  G 
x  x  x  v  v  v 


(1.2.39) 


-%  --v.  v- 


-I  'J  .  -■  v  v  J 


where  G  is  the  tide-generating  force,  2  is  the  Coriolis  factor, 

-  =  2w  sin  2  (w  =  angular  velocity  of  the  earth  rotation,  <J>  is  latitude), 

and  g  ,  g  are  Che  components  of  gravity  in  the  horizontal.  Assume  the 

x  y 

following : 

1 .  g  ,  G  <<<<  f2v  , 

X  X 

2.  g  ,  G  <<<<  2u  , 

V  V 

then , 

F  =  ftv  ,  F  =  -Qu 
x  y 

Integrating  over  the  vertical 


f 


2!v  dz  =  Qvh 


f 


-ftu  dz  =  -Quh 


(1.2.35) 


where  u  ,  v  are  vertically  averaged  velocity  components,  and  h  »  n  -  Z  . 

100.  For  a  turbulent  flow  an  eddy  viscosity  z  is  employed  in  the 
place  of  the  dynamic  viscosity  u  .  The  terms  become  using  and  for 

horizontal  and  vertical  eddy  viscosity,  respectively: 


\2  .2  \  /l 

.  o  u  .  3  u\  9  u 

Ch  -2+7T  +  Ev~7 

,3x  3y  /  3z 


(1.2.36) 


32v  2^v 

S  l  2  2  )  Ev  2 

k3x  3y  /  3z 


lnl.  The  horizontal  eddy  viscosity  terms  are  much  smaller  than  the 
vertical  eddy  viscosity  terms  and  have  been  neglected  by  some  modellers.  We 
consider  the  terras  here  in  the  following  manner.  If  we  consider  the  u 
equation  and  note  u  -  u  +  u'  ,  integrate  over  the  vertical,  and  employ 
'.eibniz  rule,  the  following  relation  is  obtained  for  the  e  terms. 
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It  is  assumed  that  all  derivative  terms  may  be  neglected,  thus 


32(u  +  u') 


32  ,,-s  .  92u  3h  3u  -  32h 

7T  (hu)  =  h  7T  +  2  37  37  +  u  7T 

3x  3x  3x 


(1.2.39) 


It  is  further  assumed  the  second  and  third  terms  are  negligible.  If  similar 
assumptions  are  made  for  the  other  terms  in  (1.2.36)  we  obtain 


,2  \ 

3  u  ,  3  u  t  , 
- r  +  - r  dz 
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h  3x2  3y 2 
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(1.2.40) 


.  .V.  *,/,  V-.V 


cv 


These  terms  are  retained  in  the  motion  equations  and  serve  to  stabilize  the 
numerical  approximations.  The  vertical  eddy  viscosity  terms  are  integrated 
over  the  vertical  as  follows: 


where 


where 


J 
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dz 
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sx 
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(1.2.42) 


(1.2.43) 


t  ,  t  =  surface  stresses 

sx  sv 

,  t.  =  bottom  stresses 

by  bx 

102.  Consider  the  bottom  stress,  x,  ,  as  follows: 

b 


T 


b 


o 


=  a  drag  coefficient 
V  =  the  fluid  velocity 


Letting  Chezv  c  =»  <^2g/C 


obtain : 


tb  =  ^  d  v 
c 


(1.2.44) 


Resolve  i,  along  the  x  and  y  directions  noting  V 


— ~> 

=  V  u“ 


— ? 

+  v“  to 


obtain : 


-  =  2S  V  u 

bx  2  f  u 
c 


Tby  =  *2  Vf  V 
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(1.2.45) 


103.  The  surface  stress  t  may  have  a  similar  form;  namelv, 

s 


p  v 

_  a  w 
Ts  '  f  2 


where 


c^  =  drag  coefficient 


o  =  air  density 
a 


v  =  wind  speed 
w 


Assuming  the  shear  stress  varies  linearly  with  depth  we  obtain 


3x 

3z 


-  +  - 
s  b 


=  ll  +  V\  _  9pwd 


as 


(1.2.46) 


where 


“wd 

s 

Integrating 


5  pressure  intensity  produced  by  the  wind 
:  distance  in  the  downwind  direction 
(1.2.46)  over  the  vertical 
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is  Che  angle  between  the  wind  direction  and  the  +  x  axis, 
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=  KV 

cos 
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sx 
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=  KV2 
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A 

sy 

w 

( 1.2.18) 
(1.2.49) 


104.  Assembling  our  results,  we  obtain  the  final  expression  for  the 
depth  integrated  motion  equations 
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Setting  S  =  1  and  expanding  the  left-hand  side  of  the  above  two  equations 
one  obtains 
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Letting  =  ot  (t  is  a  kinematic  eddy  viscosity)  and  dropping  the  bar 
n  \ 

notation  we  obtain  the  final  form  of  the  equations: 

3u  3u  3u  ^  (  3“u  ,  3“u\  u(u  +  v)1^2  ,  K  „2 
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(1.2.5 


Forcing  functions 

’.n5.  In  the  simulation  of  hurricane  induced  water  level  fluctuations, 
it  is  necessary  to  consider  the  atmospheric  pressure  anomaly  and  wind  stress 
terms  in  the  motion  equations.  The  formulations  employed  are  presented  In 
detail  in  turn  below. 


Atmospheric  Pressure  Anomaly 


If  we  define 
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=  I. 1 4  AP 


and  assume 


P  is  a  constant,  then  we  obtain: 
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3r 
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3  n 

a 

a 

a 

a 
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(1.3.5) 


If  we  employ  (1.3.5)  in  (1.3.2)  the  following  alternative  form,  which  is 
employed  in  the  current  model  is  obtained: 
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(1.3.6) 


Surface  Wind  Stress 


108.  The  hurricane  wind  loading  is  implemented  as  a  surface  stress  in 
each  horizontal  motion  equation  (Refer  to  Equations  1.2.54  and  1.2.55).  In 
the  previous  development  K  -  (c’  lo^)/ 2  and  is  dimensionless.  If  •  »  2  . 
then  K  *  P  c'  .  Several  investigators  had  developed  relationships  based 

-  3.  L 

f —  ^  — 

upon  observed  wind  data  for  c  ^  .  These  relationships  are  presented  in  turn 
below,  in  which  =  c^.  . 


Van  Dorn  Formulation 


109.  Professor  W.  G.  Van  Dorn  measured  the  wind  induced  s',  pe  u  f.e 
surface  in  an  800  foot  model-yacht  pond  [181  and  presented  the  'win; 

result . 


S  <  Kr 


.895 


,895  +  2.034(  1  -  7.2/1’  5**: 

iU 


;  o 


se  r 


se  e 


where 

C =  overwater  neutral  10  meter  drag  coefficient 
=  sustained  (10  min)  average  windspetd  at  10  m 

110.  This  formula  has  been  traditionally  used  in  hurricane  wind-stress 
.rn I i _at ions  in  numerical  storm  surge  modelling.  The  formula  was  developed 
. :  -  i  w  indspeeds 1  ess  than  15  m/sec. 
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Charnock  Formulation 


harnock  [19]  proposed  the  following  relation: 
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is  very  s 

imilar  to 

■e  'harnock  lower  Law  presented  above. 

1L.  Wu  has  developed  the  following  law  similar  to  the  Charnock  Linear 
iv  which  nav  he  used  in  hurricane  conditions  [2  L,  22]: 


I' 

'D 


«  10 


Wu  Linear  Law 


(1.4.6) 


m 

•"ere  ini  '  ,  are  as  previously  defined. 

Powell  Formulation 


6 .  Four  diagnostic  marine  boundary-layer  models  were  evaluated  for 
ipp 1  icab i  1  it v  to  the  hurricane  regime.  Model  results  also  included  the 
11  meter  level  neutral  drag  coefficients  [23].  Powell  suggests  the  following 
relation  as  an  alternative  to  the  Deacon  relationship  employed  by  Rosenthal 
'2-1  in  the  National  Hurricane  and  Experimental  Meteorological  Laboratorv 
computer  model. 


CQ  *  103  =  1.0236  +  5.366  *  10-2  U. ^  Powell  Law  (1.4.7) 

CD  *  103  =  1.1  +  4.0  x  10“2  U1Q  Deacon  Law  (1.4.8) 

where  and  1'^  are  as  previously  defined.  The  Powell  Law  has  been 

developed  based  upon  data  from  seven  hu--lcanes  and  is  the  form  which  will  be 
used  during  the  study. 


Implicit  Formulations 

116.  The  following  set  of  relationships  is  used  to  determine  the 
neutral  ten  meter  drag  coefficient. 


=  Ku1Q/ln(10/zo) 

(a) 

-  a/8  '4 

(b) 

(1.4.9) 

(c) 

where 


m 


I'*  r.  friction  velocity 
*  =  von  Karman's  constant 
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^aerodynamic  roughness 
g  ^gravity 

a  :Charnock's  constant 
and  Cp  ,  l’  are  as  previously  defined. 


Alternative  values  have  been  used  for  "■  and  a  as  shown  in  Table  III-L 
below. 


Table  III- 1 

Coefficient  Values  for  Implicit  Drag  Formulations 


nvestigator 


Surface  Wind  Stress  Evaluation 


II'.  In  order  to  compare  results  obtained  from  the  above  drag  laws  v.vr  i 
range  of  ten-meter  windspeeds  Program  DRAG  as  listed  in  Appendix  A  was 
developed  after  Prater  [25].  Results  are  shown  in  Table  1 1 1  —  2  for 
10  t  ■_  120  kts  at  10  knot  increments.  The  first  line  associated  with  a 

given  wind  speed  contains  the  drag  coefficients  predicted  bv  each  formulation. 
The  second  line  associated  with  the  given  windspeed  contains  the  ratio  of  the 
various  drag  coefficients  to  the  van  Dorn  coefficient,  which  has  long  been 
used  as  the  standard.  As  observed  the  ratios  are  all  greater  than  one  and  in 
general  increase  with  windspeed.  At  all  windspeeds,  there  are  variations  of 
over  302.  Variations  on  the  order  of  50-I00Z  may  occur  at  hurricane 
windspeeds . 


Bottom  Friction 

118.  In  Che  treatment  of  hurricane  induced  water  levels  In  hake 
Okeechobee,  it  is  necessary  to  consider  the  role  of  bottom  friction  in  Equa¬ 
tions  (1.2.54)  and  (1.2.55).  A  Chezv  0  formulation  is  employed,  which  is 
based  on  the  following  relationship  in  terms  of  Manning's  n  roughness. 

C^i^iR176 

n 

where 

C  = Chezv  C 

n*  = Manning's  n  roughness  (effective) 

R  -  hydraulic  radius  (depth) 

The  user  specifies  Manning's  n  versus  stilled  water  depth  correspondence. 
This  correspondence  is  further  modified  by  employing  the  following  relation¬ 
ship  to  account  for  the  canopy  effects  of  vegetation  near  the  shoreline. 

(n  fl  +  k  e“  d/k2>  d  ♦.  d* 

n  d  ,  d* 

w  here 

n*  -  Effective  Manning's  n  roughness 
k.  Canopy  coefficient  number  one 

r’. ,  'anopy  coefficient  number  two 

d  :  Water  depth 
d*  :  Canopy  water  depth 

Equations  (  1 .  4 .  1 )  and  (1.4.2)  are  used  to  determine  the  Chezv  C  for  each  cell 
at  each  time  step  in  the  numerical  solution.  Observe  in  this  approach  as  the 
water  depth  decreases,  n*  increases  and  C  decreases.  Since  0  appears  in 
the  denominator  in  both  motion  equations  (1.2.54)  and  (1.2.55),  the  magnitude 
of  these  friction  terms  increases. 
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onsider  the  continuity  balance,  Fquation  (l. 1.9),  the 


mronent  s  •  -  ner.tum.,  Equations  '  .  2  .  hi )  and  <1.2 .  5S)  ,  and  the  appropriate 

ressure  gr  *  riert  terns  in  Equations  f  I  .  1 .  s'*  and  (  .  3.6',  for  the  respective 

-er.f.~  e  r:  i*  :  'ns .  1  t  we  substitute  d  for  h  we  obtain  the  following  set 

•••  ;r  ;""i~i  i'  :  :.u  :  ms  written  in  subscript  notation: 


a  x 


u' 
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1  u 

XX 


a 


i  ■  ■  *  j 


w  i  mi  stressed,  etc. 

*■  >  1  n  •  a  1  rate 

r  t  ■'  1  i  s  parameter 

•  !  *.* •;  sms  i  t  v  £*>«: :  ic  ient 

O'ir’stati  elevation  due  to  atmospheric  pressure  anora v 

i 

water  sur:  ace  elevation  with  respect  to  tire  lake  datum 
tto-  elevation  with  respect  to  the  lake  ■latum 

:  •’  -  pi.it  .  -r  !  ‘  ,  the  right  hand  side  has  been  set  to  the  rain  tall 

-  •  e  f  equation  set.  dote  also  that  in  Fquation  i  1  .  i .  hi  fa)  and 


r 
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j 
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ud 


V  cos  d 

w 


F 

v 


V4 

ud  w 


s  in 


a 


where  J  is  the  densitv  of  air. 
a 

Voting  that  0  corresponds  to  either  u  or  v  or  r.  we  rewrite  Equation 


r  1  .  5  .  Ft  a  in  transformed  form  as  follows: 


+ 


I  /  -0 


(1.5.9) 


t 


fu  +  g  -  “  +  (u*-  +  v  ) 

u2  C“d 


(b) 


(c) 


(  1.5.9')  is  the  object  of  numerical  approximation.  In  the  initial 
a  linear  form  of  the  equations  is  considered  by  setting  all 
eras  to  zero.  In  addition,  ^  ,  and  d  are  all  considered 


Short  Wave  Equations 

In  order  to  aid  in  the  development,  let  us  define  the  following 
ers . 

Wave  height  measured  from  trough  to  crest 

Significant  wave  height,  defined  as  the  average  height  of  the 
highest  one-third  of  the  waves  recorded  in  a  specific  time 
period 

Wave  length  from  crest  to  following  crest 

Wave  period  (time)  from  crest  to  crest 

Average  depth  of  water  in  the  wave-generating  area 

Fetch  or  wave-generating  distance 

Windspeed  over  the  fetch 


t  Duration  of  the  wind  over  the  fetch 


MWL  Mean  water  level 

VTL  Vind  induced  level 

T  Significant  period,  defined  as  the  dominant  period  in  the  groups 
of  highest  waves  recorded  in  a  specific  time  period 

g  Gravity 

Semi-empirical  equation  development 

122.  Traditional  approaches  have  involved  the  development  of  semi- 
empirical  equations  for  and  Tg  in  terms  U  ,  F  ,  t  ,  and  d  .  In 

deep  water,  the  so-called  SMB  Method  after  Swerdrup,  Munk,  and  Bretschneider 
[26]  is  embodied  in  the  following  three  equations  presented  in  the  Shore 
Protection  Manual  (SPM)  [27]. 


0.283  tanh 


0.0125 


Si 

..2 


0.421 


(2.1.1) 


1.20  tanh 


(2.1.2) 


and  , 


where 


exp 


and 


•  x  • 

x  •  =  e 
In  =  log^ 

K  =■  6.5882 
A  =  n.0161 
3  -  0.3692 
C  -  2.2024 

D  -  0.8798 
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Given  U  and  t  ,  solve  Equation  (2.1.3)  for  F  and-s6t  F  =  F  .  Then  compute 

m 

F  eff  =  min  (F  ,  F  )  ,  where  F  is  given.  If  F  eff  =  F  ,  the  wave  growth  is 
m 

duration  limited.  If  F  eff  =  F  ,  then  wave  growth^is  fetch  limited. 

m 

Finally,  erolov  1.1.  L)  and  (2.1.1)  with  l'  and  F  eff  to  compute  H  and 

s 

123.  In  the  above  method,  a  constant  windspeed  is  assumed  to  act  for  a 
duration  t  over  a  fjetch  F  to  generate  the  deepwater  waves.  The  wave 
growth  can  be  limited  by  fetch  length  or  by  wind  duration. 

124.  More  recently,  work  has  focused  on  the  development  of  directional 
wave  spectral  models,  in  which  an  energy  distribution  is  determined  versus 
frequency  for  each  direction.  Significant  wave  height  and  period  are  directlv 
determined  for  each  direction  based  upon  characteristics  of  the  energy  versus 
frequency  distribution. 

125.  The  directional  wave  spectrum  approach  is  adjudged  to  be  outside 
the  limits  of  project  resources.  The  SMB  method  has  also  been  extended  to 
shallow  water  as  described  in  the  following  two  equations  presented  in  the  SPM 
[2],  in  which  a  friction  factor  f^  =  .01  and  percolation  effects  have  been 
considered . 


Mote  H  and  T  are  determined  in  terms  of  F  ,  U  ,  d  ,  and  f,  *  .01  . 
s  s  f 

126.  Recently  Equation  (2.1.1)  through  (2.1.5)  have  been  updated  based 
upon  wave  spectral  model  studies  (CETN-I-7  Equations  1-3  [29]  and  CETN-I-6 
Equations  4-5  [28]).  Since  conditions  in  Lake  Okeechobee,  relate  to  shallow 
water  wave  generation,  the  updated  shallow  water  relations  are  presented  onlv 
as  follows. 


12'.  Note  in  the  above  equations  the  windspeed  U  has  been  replaced  by 

an  adjusted  windspeed,  U  .  The  following  assumptions  are  made  in  deter- 

A 

mining  U  as  outlined  in  CETN-I-5  f30] : 

t\ 

a.  The  windfields  are  well  organized  and  can  be  adequately 
described  by  the  use  of  an  average  windspeed  and  direction 
over  the  entire  fetch. 

b.  The  windspe  d  should  be  corrected  to  the  10  meter  level. 

c.  The  windspeed  should  be  representative  of  the  average  wind- 
speed  measured  over  the  fetch.  That  is,  for  long  fetches  with 
more  than  one  wind  station  the  windspeed  is  averaged  for  the 
stations.  The  windspeed  is  further  adjusted  for  the  minimum 
duration;  the  average  windspeed  for  the  duration  is  determined 
from  the  wave  forecasting  curves.  This  process  may  need  to  be 
repeated  to  obtain  the  final  adjusted  value  of,  U  . 

A 

d.  The  windspeed  must  be  adjusted  to  account  for  the  non-constant 
coefficient  of  drag. 

e.  «'hen  the  fetch  length  is  10  miles  (16  km)  or  less,  the  wind 
has  not  fully  adjusted  to  the  frictional  characteristics  of 
the  waves.  In  such  cases,  the  overwater  windspeed  will  be 
estimated  to  be  110  percent  of  the  overland  windspeed,  U  • 

Thermal  effect  on  stability  of  the  air  in  this  case  is  not 
applicable . 

f.  When  the  fetch  length  is  greater  than  10  miles  (16  km), 
thermal  stability  effects  must  be  included  in  the  windspeed 
transformation . 
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to  be  made  at  hourly  intervals.  Equation  i  2  .  1  .  ?  '  wi L 
wind  duration.  t  .  For  compatibility  with  Equation  _  .  1  .  ~ 
(.33,  1.5)  hours  as  indicated  in  Figure  I  in  i  Since  -  • 
usuallv  less  than  30  feet  in  bake  Okeechobee  even  under  st  r- 
growth  may  be  limited  by  depth. 

CW-lb?  Curves 

L29.  A  second  approach  utilizes  curves  developed  Jurir 
sive  Project  CW-167  [  3  1  j  on  bake  Okeechobee.  Figures  .;o  and 
represented  by  the  following  relation: 


'e  y  . 

1  .  i  "  t* 

1  v  e  4  * 

i  e  r  t 


-ei  *■ 


\ r 


: '  e 


H 


s  -  g  man 


.003563252 


.  >448 1  5  352  .  h Q 0 0 Q  l  l  " * 

Si  \  t  . 18438944“!*  /  Si 


UA 


where  ,  F  ,  d  ,  and  g  are  as  defined  previously 

Figure  lb  in  [31 ]  is  represented  as  follows: 


f  .9056838' 

T  =  b.  2620  18: 14  —  (  ££ 

e  W 


where  T  ,  U  ,  d  ,  and  g  are  as  previously  defined, 


'•'hve  runup 

130.  Wave  runup  has  been  measured  in  the  laboratory  using  physical 

models  with  respect  to  mean  water  level  and  as  a  result,  wave  set-up  is 

implicitly  included.  No  separate  analysis  is  developed  here  to  distinctly 

separate  out  wave  set-up  from  the  total  runup.  Consider  Figure  1 1 1 —  1  taken 

from  Stoa  [32)  to  define  the  variables  applicable  to  runup.  The  beach  slope, 

r  ,  is  used  in  conjunction  with  Figure  7-4  in  the  SPM  [27],  to  determine  the 

maximum  wave  height  which  may  be  expected  to  break  on  the  beach  prior  to 

reaching  the  structure.  If  Hg  as  previously  computed  is  greater  than  the 

maximum  breaker  height,  H  ,  ,  H  is  set  equal  to  H  ,  .  The  modified 

mb  s  mb 
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is  computed  using  Figure  7-5 


deepwater  waveheight,  H'  ,  associated  with  H  , 

o  mb 

in  the  SPM  [2 7 ] . 

131.  Procedures  developed  by  Stoa  [32]  for  a  horizontal  nearshore 
bottom  are  employed  as  shown  in  Figure  III-2.  Equation  8  in  Stoa  [32]  as 
referenced  in  this  figure  is  as  follows: 


H 1  \ q—  ^ 

ttt  -  (cot  6)‘1,04  (4.23)  102(q_1)  ~f  for  cot  0  >  2  (2.3.1) 

o  VgT  / 


where  variables  are  as  indicated  in  Figure  III-l.  In  order  to  determine  q  , 
Figure  4  in  Stoa  [32]  is  utilized  for  cot  9  =  6  or  8  ,  which  represent  the 
range  of  levee  slope  conditions  at  Lake  Okeechobee.  All  computed  values  of 
runup  are  multiplied  by  values  in  Figure  50  in  Stoa  [ 32 J  to  account  for 
physical  model  scale  effects. 

Wave  overtopping 

132.  Weggel  [33]  based  upon  an  analysis  of  physical  model  results  has 
presented  the  following  equation  to  predict  wave  overtopping: 
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overtopping  rate  (L“/T)  per  unit  length  of  levee 
gravity 

2 

overtopping  rate  (L  /T)  per  unit  length  for  a  levee  with  crest 
elevation  at  the  stilled  water  level  (empirical) 
deep  water  significant  wave  height 

runup  height  measured  vertically  from  the  stilled  water  level, 
e.g.,  the  height  to  which  the  water  would  runup  if  the  levee  were 
high  enough  to  preclude  overtopping 

height  of  the  levee  above  the  bottom 

water  depth  at  the  levee  toe 

empirical  coefficient 
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Figure  1 1 1  —  2 .  Lake  Okeechobee  Runup  Mechanics  (after  Stoa  [32]  Figure  6) 

L33.  Note  ct  and  Q*  are  empirically  determined  coefficients  which 

o 

depend  on  incident  wave  characteristics  and  levee  geometry.  Approximate 
values  of  a  and  Q*  for  various  slopes  and  structure  types  are  given  in 
Figures  7-24  through  7-32  in  the  SPM  [27]  as  function  of  wave  steepness, 


x*r **“ « 


twgT  ,  and  relative  height,  dg/H^  •  Weggel  [33]  has  suggested  that  the 
following  two  relations  be  employed  in  cases  where  more  exact  values  are  not 
available.  His  procedure  uses  theoretical  results  for  wave  overtopping  on 
smooth  slopes,  assumes  a  1:10  beach  slope  as  do  Figures  7-24  through  7-32,  and 
yields  conservative  overtopping  rates. 


a  =*  .06  -  .0143  In  (sin  9) 


2  ,  2ird 


(2.4.2) 


rp  tanh  (  __s  ) 


where 

6  £  levee  slope 
L  =  wave  length 
Hg  =  significant  wave  height 
e  £  parameter  based  upon  wave  slope 
T  =  wave  period 

a  ,  Q*  ,  ,  dg  ,  and  q  are  as  previously  defined. 

134.  For  a  sinusoidal  wave,  e  *  l/2n  ,  and  assumes  a  maximum  value 
relative  to  a  coidal  wave  shape.  The  above  relations  with  z  -  1/2t  are  used 
in  the  numerical  model.  Equation  (2.4.1)  does  not  need  to  be  corrected  for 
scale  effects  if  the  scale-corrected  runup,  R  ,  is  used  in  the  equation. 
Onshore  wind  effects  increase  the  overtopping  rate  at  a  given  levee  section. 

As  a  guide,  the  following  multiplication  factor  is  presented  in  the  SPM  [27]: 


h  -  d 


V  R 


(2.4.3) 


where 


k'  £  multiplication  factor  due  to  onshore  winds 

w^  £  onshore  windspeed  coefficient 

h  ,  d  ,  R  and  9  are  as  previously  defined, 
s 


i, 


The  following  relations  are  suggested  for  : 


w  —  •  5  V  /  3  0 
f  ws 

w  =  .5  +  1.5  (V  -  30) / 30 
f  ws 

w,  =  2 


0  ^  V  ^  30  mph 
ws  r 

30  <  V  ^  60  mph  (2.4.4) 

ws 

V  >60  mph 
ws  r 


where  V  is  the  onshore  wind  in  miles  per  hour.  This  multiplication  factor 
is  used  in  the  numerical  model  in  order  to  make  the  most  conservative 
overtopping  rate  estimate. 

135.  Equation  (2.4.1)  assumes  that  the  levee  crown  elevation  is  greater 

than  the  stilled  water  elevation;  i.e.,  0  <  (h  -  d  )/R  <  1  .  Let  us  con- 

s 

sider,  the  two  cases  not  included.  If  (h  -  dg)/R  ^  1  ,  then  Q  *  0  ;  i.e.,  no 
overtopping  occurs.  If  (h  -  dg)/R  2  0  ,  the  situation  is  as  indicated  in 
Figure  III-3. 


Figure  III-3.  Overtopping  Case  for  (h  -dg)/R  ^  0 

136.  Consider  a  surface  wave  represented  by  W(x)  -  Hg/2  sin  21’x/L  , 

where  x  is  a  distance  along  the  wave.  The  elevation  of  the  water  surface  is 

then  we(x)  =  w(x)  +  d  .  If  we  assume,  d  -  h  ^  H  .'2  ,  then  we  seek  to  find 

3  s  s 

xq  ,  such  that  we(xQ)  *  d  -  h  .  The  following  relation  may  be  used. 


(2.4.5) 


where  all  quantities  are  as  previously  defined.  To  obtain  the  volume  of  water 
passing  over  the  levee,  the  following  approach  is  utilized: 


PART  IV:  HYDRODYNAMIC  MODEL:  NUMERICAL  IMPLEMENTATION 


Long  Wave  Equations 

138.  The  hydrodynamic  model  consists  of  a  long  and  short  wave  compon¬ 
ent.  The  long-wave  equations  as  previously  developed  are  solved  on  a  space- 
staggered  grid  employing  an  alternating  direction  implicit  (ADI)  scheme.  The 
details  of  this  numerical  implementation  are  presented  initially.  In  the 
second  section,  the  numerical  approach  used  to  implement  the  previously 
developed  short  wave  equations  is  presented.  The  essential  ingredient  in  the 
numerical  implementation  is  the  methodology  for  determining  fetch  in  curved 
hurricane  windfields.  The  numerical  implementations  are  presented  in  turn 
below. 

Development  of  the  stabilizing  correc- 
tion  scheme  for  the  long-wave  equations 

139.  In  order  to  develop  the  finite  difference  approximations,  first 
consider  the  following  reduced  form  of  equation  set  (1.5.9)  in  which  all 
underlined  terms  are  set  to  zero. 

u  +  ^—  n  -  0  (a) 

t  ^ 

u  +  n  »0  (b)  (1.1.1) 

t  u2  a2 

n  +  —  u  +  —  v  =0  (c) 

t  uL  u,  u2 

The  above  equations  essentially  ignore  surface  and  bottom  stress,  Coriolis 
force,  convective  acceleration,  lateral  diffusion,  rainfall,  and  atmospheric 
pressure  anomaly.  Observe  further  that  the  depth  is  also  considered  constant. 
Stretching  functions  and  u2  are  also  constant. 

140.  Equation  set  (1.1.1)  is  written  in  vector  form  as  follows: 


u  +  Au  +  Bu  =  0  (1.1.2) 
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In  what  follows  the  transformed  position  and  time  coordinate  (a^,  t)  is 

represented  on  the  space  staggered  grid  shown  in  Figure  IV- 1  by  (nAc^,  mAct^, 
kAt)  .  Note  all  grid  dimensions  are  uniform  in  the  transformed  plane  grid. 

The  water  surface  elevation,  V  ,  field  is  computed  over  the  cell  centered 
field.  The  y  -  a  direction  velocity  component  field  is  computed  over  a 
field  defined  by  staggering  the  cell  centered  field  to  right  half  a  grid  cell. 
Similarly,  the  x  -  direction  velocity  component  field  is  computed  over  a 
field  defined  by  staggering  the  cell  centered  field  to  the  bottom  half  a  grid 
cell . 

141.  As  a  prelude  to  the  actual  formulation  of  the  difference  equa¬ 
tions,  let  us  consider  the  following  definitions: 


=  water  surface  elevation  at  (n,  m)  at  time  level  k  +  I 

=  -  velocity  component  at  (n,  m+l/2)  at  time  level  k  +  I 

=  -  velocity  component  at  (n+1/2,  m)  at  time  level  k  +  I 

e  (-1,  0,  L) 

Observe  in  Figure  I V— 1 ,  the  still  water  depth  field  h  ,  is  specified 
on  the  cell-centered  grid.  The  depth  convention  employed  is  shown  in  Fig- 

k 

ure  IV-2  for  conditions  at  time  level  k  .  It  is  assumed  h  >>  n  and 

n ,  m  n ,  m 

d  '  -  -h  a  constant  for  all  time  levels  k  . 
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142.  In  developing  Che  general  approach,  consider  Che  following  chree 


Cime  level  finice  difference  scheme: 


u  -  u  ,  <Sx  j  „  ,  k+1  ,  k—  1 N  /  .  6y  k+1  ,  k-1. 

- -  — -  +  —  ■  A(u  +  u  )j  +  -4  I  B(u  +  u  ) 


=  0  (1.1.3) 


where 


n , m+ 1/2 


n+1/ 2 ,m 
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n  ,m+I _ n  ,m 
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n+l,m  n,m 
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5xu^  = 


un,m-l/2  un,m-l/i 
UiAa! 


Vn+l/2,m  Vn-l/2,m 
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Observe  from  Equacion  (1.1.2)  ChaC  in  macrix  A  associaCed  wich  u  Che 

ai 

chird  column  concains  all  zeros.  Therefore,  no  encry  is  required  in  Che  chird 

row  for  5x  .  Similarly,  in  macrix  B  associaCed  wich  u  ,  Che  second  column 

a  2 

concains  all  zeros.  Thus,  no  encry  is  required  in  5y  for  Che  second  row. 
143.  If  we  define  Che  operaCors  a  *  ACdxA  and  X  =  4c5yB  Co 

k+ 1  x  y 

operace  u  ~  from  lefC  Co  righc,  we  may  rewriCe  (1.1.3)  in  che  following 
form: 


(I  +  X  +  A)u 
x  y 


(I-X  -X)u 
x  y 


(1.1.4) 


k+1  k-l 

144.  Adding  X  X  (u  -  u  )  Co  Che  lefc  hand  side  of  che  above 

x  y  k+i 

equation  we  obcain  a  factorizaCion.  The  addicion  of  the  term  X  X  (u 

k-i  x  y 

u  )  does  not  make  the  difference  approximation  inconsistent  and  the 

2  2  2 

truncation  error  remain'-  0[At  ,  Ax  ,  Ay  ]  for  =  1  . 


145. 


We  have  thus  achieved  an  approximation  of  the  following  form: 
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v.  -■ 


( I  +  ’>  )(I  +  )uk+L  =  (I-x)(I->.  )uk  1 

XV  x  y 


(1.1.5) 


If  an  intermediate  solution  level  u*  is  introduced,  there  are  many  possible 
schemes,  which  may  be  used  to  split  the  above  equation.  The  present  formula¬ 
tion  employs  a  stabilizing  correction  scheme  as  given  below. 


(I  +  >  )u*  =  (I  -  >  -  2>  )uk_1  (a) 

X  XV 


/T  ,  .  k+l  A  ,  k-1 

( I  +  \  )  u  =  u  *  +  X  u 

y  y 


(1.1.6) 


Equation  (1.1.6) (a)  is  treated  in  the  -  x  sweep  ,  while  equation 

(1. 1.6)  (b)  is  handled  in  the  do  _  y  sweep  .  Hence,  the  name  multi- 

operational  (two  sweeps)  alternating  direction  (first  -  x  then  -  v) 

implicit  finite  difference  method  of  stabilizing  corrections. 

146.  In  order  to  form  the  foundation  for  the  approximation  of  the 
complete  nonlinear  hydrodynamic  equations,  let  us  first  expand  equation  set 

(1.1.6) 


x  sweep  (Equation  1.1.6  (a)) 
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147.  In  further  expanding  the  above  equation,  the  difference  operators 

;x  and  5v  are  applied  to  the  matrix  vector  products.  The  A  and  A  are 

xv 

performed  on  u  from  right  to  left.  Therefore,  we  obtain  for  the  x-sweep: 
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j,  -  y  sweep  (Equation  1.1.6(b)) 
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(1.1.8) 


(1.1.9) 


Expanding  the  above  equation  for  each  component,  we  obtain: 
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If  we  now  substitute  (1.1.8)(c)  into  (1. 1. 10) (c)  and  substitute  (1.1.10) (b) 
into  (1.1.8) (a)  and  (b)  we  obtain  the  following  forms  for  the  two  sweeps. 


x-sweep 
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(1.1.11) 


y-sweep 

k+1  ,  itd  k+1 
+  - - -  v 


k+1 
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=:  r* 


Atd  k-L 
+  - : -  v 


k- 1 


U-,Aa  n+l/2,m  n-l/2,m 
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k+l  Atg  /  k+1  k+1  \  k-1 

n+l/2,m  ,  n+l,m  n,m/  n+l/2,m 

*"  \  / 


(a) 


Cl. 1. 12) 


_  _AIS_  /  nk_l  -  nk*L 
i  n+l  »m  n.ni, 


(b) 


148.  If  In  Equation  sets  (1.1.11)  and  (1.1.12),  all  terms  on  the 

right-hand  sides  are  brought  over  to  the  left-hand  sides  and  the  resulting 

eauations  are  divided  by  2it  ,  the  finite  difference  approximations  are 

obtained  in  a  more  standard  form.  Equation  (l.l.ll)(a)  represents  the 

approximation  to  the  continuity  equation  (Equation  (1.1.1) (c  )),  Equation 

(l.!.ll)(b)  represents  the  approximation  to  the  a[  -  x  motion  equation. 

(Eauation  (l.l.l)(a)),  and  Equation  (1.1.12Xb)  represents  the  approximation  to 

the  -  v  motion  equation  (Equation  (l.l.O(b)).  Equation  (1.1.12)(a)  is  a 

k+1 

tabilizing  correction  to  obtain  n  from  "*  and  involves  only 

n  ,m  n  ,m 

■-‘t-c opponent  velocity  terms. 

149.  We  are  now  in  position  to  obtain  approximations  for  the 

nonlinear  hydrodynamic  set  (Equation  set  (1.5.9)).  We  do  so  in  a  two  step 
process.  In  step  one.  we  relax  the  assumption  -  j  =  =  1  and  that 

d„  =  -h  _  .  Equations  (1.1.11)  and  (1.1.12)  now  become: 


150.  Observe  in  Equations  (1.1.11)  and  (1.1.12),  all  stretching  func¬ 
tions  and  u7  are  evaluated  at  the  center  of  the  appropriate  differences 

shown  in  their  respective  parentheses.  The  depth  is  evaluated  at  time  level 

k  at  the  location  of  the  corresponding  velocity  as  follows: 

-  /2 

-  hn,„)  /2 

151.  In  step  two,  the  previously  deleted  (underlined)  terms  in  Equation 
set  (1.5.9)  are  considered.  In  the  continuity  equation,  the  rainfall  term  is 

k 

introduced  in  Equation  (1.1. 11) (a)  by  adding,  2itR  ,  to  the  right  hand 
,  m,n 

side,  where  K  is  the  rainfall  rate  at  (n,m)  at  time  level  k  .  Note 
m,n 

Equation  (1.1. 12) (a)  is  not  changed.  For  the  -  x  motion  equation,  the 
convective  acceleration  (l) ,  Coriolis  @  ,  atmospheric  pressure  anomaly  (3) , 
bottom  friction  ©,  lateral  momentum  diffusion  0,  and  wind  stress  ©  terms 
are  approximated  in  the  following  manner  to  obtain: 


n±l/2,m 


n,m±l/2 


'  k  u  .  k 

n  -  h  ,  +  n 

ni  1  ,m  r.t  1  ,m/  \  n ,m 


^n,m±l  ^n,m±l^  +  ^nn,i 

m 


The  approximation  of  each  term  is  discussed  briefly  below: 


I 

| 

t 


I 


(p  Convective  acceleration  term 

All  differences  are  centered  in  space  and  in  time  (at  time  level  k)  in 
order  to  maintain  stability.  Since  the  -  x  motion  equation  is  evaluated 
at  (n,m+l/2)  and  -  x  velocity  components  are  available  only  at  (n+l/2,m), 
it  is  necessary  to  define 

=k  /  k  ,  k  ,  k  ,  k 

Vn,m+1  '2  \Vn+L/2,m  Vn-l/2,m  Vn+l/2,m+l  Vn-l/2, 

(2)  Coriolis  term 

s{( 

The  approximation  is  centered  at  time  level  k  with  v^  m+\/2  as 
previously  defined.  Note  f  is  equivalent  to  Q  as  defined  under 
Equation  (1.2.34)  in  PART  III. 

CD  Atmospheric  pressure  anomaly  term 

The  approximation  is  centered  in  time  and  in  space  with  defined 

above  Equation  (1.3.5)  in  PART  III. 

(p  Bottom  friction  term 

Note  a  Chezy  C  approach  is  used  as  given  by  Equation  (1.4.1)  in  PART  III 

and  is  a  time  dependent  quantity  based  upon  water  depth.  Since  water  depth  is 

available  at  (n,m)  and  the  a  -  x  velocity  component  is  calculated  at 

/  k  1  k  \  / 

(n,m+l/2),c  , ,  =  I c  ,,+c  /2  .  It  is  also  necessary  to  empl ov 

^  n,m+l/2  y  n,m+l  n,m// 

v  as  defined  previously  with  k-1  substituted  for  k  .  The  i  -  x 

n  f  I  /  c  i 

velocity  component  is  taken  at  time  level  k+1  ,  while  the  magnitude  of  water 
velocity  is  computed  using  quantities  of  time  level  k-1  .  This  form  of  the 
friction  term  has  been  used  by  several  investigators  and  is  generally  believed 
to  be  the  most  stable  form  of  approximation  for  this  term. 

(5)  Lateral  diffusion  term 

Observe  all  differences  are  centered  in  time  at  level  k  and  in  space. 
All  stretching  functions  are  evaluated  at  m+1/2  .  All  differences  in 
stretching  functions  are  centered  at  m+1/2  . 

©  Wind  stress  term 

This  term  is  evaluated  at  time  level  k  and  is  defined  as  follows 
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1c  a  /  k  k\  k  /kk|  k  /  1 

F, -  lco”K  }  vn ,m  +  \cowx  /  "„,■*!/££  „„ 

1  w  '  '  n  ,m  v  !  n , m+ 1  '  n , m+ 1 / 2 


(1.1.15) 


where 


P  £  densitv  of  air 
a 


P  £  density  of  water 
w 

c  £  drag  coefficient  at  (n,m)  at  time  level  k 
n  ,m 

wx  £  a,  -  x  wind  velocity  component  at  (n,m)  at  time  level  k 
n,m  l 

v  £  wind  magnitude  at  (n,m)  at  time  level  k 
n  ,m 

152.  Let  us  now  consider,  the  -  y  motion  equation;  e.g.,  Equation 

(1.5.9)(c)  .  In  analogy  to  the  a  -  x  motion  equation,  the  convective 
acceleration  0 ,  Coriolis  Q)  ,  atmospheric  pressure  anomaly  0  ,  bottom 
friction  ©  ,  lateral  momentum  diffusion  © ,  and  wind  stress  ©  terms  are 
approximated  to  obtain: 


KM 


J\ 


-It  \ 


n+l/2,m  2u  ^ Aa  ^  n+i  /  2  ,m  n+l  /  2  ,m+l 


k.  Ik  k  k  =k 

Vn+l/2,m-l;  +  2  ^  A x Vn+l/2,m  Vn+3/2,m  Vn-l/2,mi+  ^Un+l/2,m 

“  **  s  / 


-  2  (-  ) 
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k+1  k-! 

+  r 
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n+1 12 ,m 
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- -Ci/2,.) +  rr^  (y  -\k) 

)  2“lMl  D  ‘Vh/2  V  1 
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n+1 / 2 ,m- 


|  (U2Aa2)  (V 


k  k  «  k 
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-  F  -  0 
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(1.1.16) 


Analogous  comments  regarding  these  approximations  hold  with  u,a.  being 

k  1  k 

replaced  by  v,  etc.  In  Equation  (1.1.15),  wx  is  replaced  by  wy  and 


n , m+ 1 / 2  n+ 1 ! 2 , m 


( 


l 


153.  Assembling  all  results,  the  final  sweep  equations  may  be  written 
in  the  following  manner. 


x  sweep 


"Vl/2  Vm-1/2  +  n£,m  +  Vl/2  Un,m+l/2  =  Am  Continuity  (a) 
-  ^+1 

-am"n ,m  +  Vl/2  uu,m+l/2  +  Vi  nn,m+l  =  Vl/2  Momentum  W 


(1.1.17) 


where 
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154.  Alternatively,  the  sweep  equations  may  be  expressed  in  operator 
form  as  shown  in  Table  IV— 1 . 

Tridiagonal  matrix  solution 

155.  It  is  instructive  to  study  the  computational  structure  embodied  in 

Equation  set  (1.1.17).  Consider  the  computational  stencils  for  the  continuity 
and  momentum  equations  as  shown  in  Figures  IV-3  and  IV-4,  respectively. 

The  locations  of  the  pertinent  hydrodynamic  variables  are  as  shown  in  Fig¬ 
ure  IV-1.  The  computational  stencils  show  the  dependence  of  the  computations 
on  the  surrounding  hydrodynamic  variables  at  different  levels  in  time.  Ob¬ 
serve  in  both  Figures  IV-3  and  4,  that  for  a  given  m  ,  the  variables  at  the 
most  advanced  time  level  all  lie  on  a  straight  line  of  constant  n  .  Since  at 

the  most  advanced  time  level,  the  variable  of  interest  (either  r*  or 

'  n  ,m 


n,m+l/ 


is  a  function  of  its  nearest  neighbors  along  constant  n  (  either 


u  ,  or 
n,m  1/2 


and  n 


n,m+l) 


,  the  difference  schemes  are  implicit, 


156.  The  solution  procedure  for  the  x  -  sweep  is  for  each  n  , 
determine  the  number  of  computational  segments  and  their  beginning  and  end 
locations  in  the  grid.  The  number  of  computational  segments  and  their  loca¬ 
tions  beginning  and  end  are  a  function  of  the  grid  system  and  the  boundary 

conditions.  In  order  to  illustrate  the  procedures,  consider  r.  -  N  and  one 

k+1 

computational  segment  defined  by  specifying  u^  and  nN  ME+1  ' 

matrix  equation  assumes  the  following  format. 


&  4-  tv  1  L  _  .  fV  T  1 

nNfMS  +  aMS+I/2UN,MS+l  /2  _  MS  ^S- 1 /2UN ,MS- 1 /2 


*  +  ~  k+1  ,  _  „ 

aMSrN,MS  aMS+l/2UN,MS+l/2  ^S+l'VMS+l  MS+1, 


k+1  .  ,  k+1 

'aJ-l/2UN,J-l/2  +  nN,J  aJ+l/2Un, J+l/2  AJ 


(1.2.1) 


(1.2.2) 


(1.2.3) 


v'vl*  //  /:■> 


*  m  +  s  -  * 

.  ^  -  n  - 

*  nJ* 


• 1 


I 


Table  IV- 1 


Finite  Difference  Approximation  In  Operator  Form 
of  the  Complete  Hydrodynamic  Equations  Set 


Fn,m+(I/2)  Fn,m-(l/2) 


2^‘n,m'  ^n+(l/2),m  ^n-(l/2),m 


; ,MF  )  =  F  -  F  =  5a  5a  (F  ) 

-I  n ,  m  n , tn+l  n ,  m- 1  1  1  n ,  n 


i,,a0(F  )  -  F  ,  -  F  =  5a  (F  ) 

2  2  n,m  n+I,m  n-l,m  2  2  n,m 


11F  =F  +F  -2F  =  6  a  5  a  F 

1  1  n,m  n,m+l  n,m-l  a,m  2121  n,m 


■  a  „  a  F  =  F  +  F  -  2F  =  6  a  5  a^F 

_  _  n,m  n+l,m  n-l,m  n,m  2222  n,i 


a,  F  ,  +  F  a  F  +F 

1  n  ,m+l _ n  ,m  2  n+l,m  n,m 


n  ,m 


n  ,m 
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(Continued) 
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"  ,m 

^  • 

• 

•  TIME  LEVEL  K 

■  TIME  LEVEL  •  OR  K  +  1 

*  TIME  LEVEL  K  -  1 


Figure  IV-3.  Computation  stencil  for  the  x  -  a.  sweep  continuity  equation 


•  time  level  k 

■  TIME  LEVEL  •  OR  K*1 
A  TIME  LEVEL  K  -  1 


Figure  IV-4.  Computational  stencil  for  the  x  -  i. 


sweep  momentum  equation 


‘aJrN,J  +  aJ+l/2Un,  J+l/2  +  aJ+lnN,J+l  =  B  J+l  /  2 


(1.2.4) 


k+1  +  *  +  k+1  -4 

'aME-l/2UN,ME-l/2  nN,ME  ®ME+1  /  2UN  ,ME+1  /  2  =  Ne 


(1.2.5) 


—  K  '  L 

~aMEn^,ME  +  aME+l/2UN,ME+l/2  =  BME -1/2  "  aME+lnN,ME+l  0-2.6) 

Observe  the  tridiagonal  structure  of  the  above  matrix  for  n  =  N  .  Let  us 
seek  the  following  general  structure 


k+1 

n*  =  -P  u*  1  .  +  Q 

n,m  m  n,m+l/2  m 


u  ,  .  =  -R  n*  +  S 

n,m+I/2  m  n,m+l  m 


m  =  MS . ME 

n  “  N 


0.2.7) 


157.  In  order  to  determine  P  ,  Q  ,  R  ,  and  S  ,  first  rewrite 

m  m  m  m 


Equation  (1.2.1) 


"k  as 

XmS  =  _aMS+l/ 2UN,MS+I/2 


k+1  r  k+1 

M.MS+I/2  +  l  MS  aMS— 1 / 2UN , MS— 1/2 


(1.2.8) 


r<>  i 

PMS  =  ^MS+l / 2  ’  QMS  ^lS  +  aMS-l/2uN, MS-1/2 


158.  Let  us  next  substitute  (1.2.8)  into  (1.2.2)  to  obtain: 


k+1  \  -  k+1 

-aMS^  MSUN,MS+l/2  %S J  aMS+l/2liN, MS+1/2 


(1.2.9) 


+  ^S+^N.MS+l  "  BMS+l/2 


159.  Rearranging  the  above  equation,  we  obtain: 
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_aMS+lnN,MS+L  +  BMS+l/2  +  a.MSQMS) 


N.MS+1/2 


aMS-M/2  + 


S,s*r— ^ - 

^S+l  /  2  +  aMSPMS 


(1.2.10) 


BMS+l/2  +  *MSQMS 


MS 


^S+l/Z  +  aMSPMS 


160.  In  order  Co  obtain,  the  general  form  for  P  ,  Q  consider  the 

m  m 

next  equation;  e.g., 


K  »  1 

'aMS+L/2(“RMSnN,MS+l  +  SMS^  +  nN,MS+l  +  ^5+3/2^, MS+3/2  =  ^SfS+l 


Rearranging  this  equation,  we  obtain: 
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MS+1  1  + 


^+3/2 
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(1.2.11) 


Ns+l  +  aMS+l/2S> 
1  +  aMS+l/2RMS 


'MS+1  1  + 


We  therefore  hypothesize  the  following  general  formulas: 
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3  a  - 

m  .  *n 
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(1.2.12) 


161.  Observe  the  above  formulas  may  need  to  be  modified  for  P  and 

n 

Q  for  m  =  MS  ,  and  for  Rm  and  Sm,n  =  ME.  Returning  to  the  and  Qm  ' ir 

x .  ^  k+1  .  . . .  ,  T  i  . 

m  =  .'lb,  no ce  LS  sPecitied.  Let. .us  rewrite  uhis  in  the  following  :  ''r 

k+1  k+1 

UM ,MS- 1/2  ”  'StS-fN.MS  +  SMS-1  '*  ^S-l  =  °  ’  SMS-1  =  UN, MS-1/2  d-2.13) 

If  we  define  and  SMS-I  as  aBove  c^en  c^e  general  equation  holds  for 

m  =  MS  as  well.  Let  us  now  consider,  the  last  equation  (Equation  1.2.6)  and 
assume  the  general  solution  form  for  : 


,_p  k+1  \  -  k+l 

'^El  MEUN,ME+l/2  gME  )  aME+l/2UN,ME+l/2 


?ME+l/2  aME+lnN,ME+l  d-2.14) 


N.ME+1/2 


‘^E+rN.ME+l  +  BME+l/2  +  ^E^E 


^+1/2  +  aMEPME 


Therefore 


^lE+I  /  2  +  aMEPME 


_  BME+I/2  +  Ne^E 
5me  "  - 

aME+I/2  ^E  ME 


We  note  that  the  standard  formulas  hold  for 


and  S._  .  We  therefore  have 
ME 


developed  the  procedure  outlined  in  Table  IV-1  for  column  n  in  the  grid. 

Step  2  is  the  forward  elimination  step  and  Step  3  is  the  backward  substitution 
step  of  the  Thomas  algorithm. 

162.  The  treatment  of  other  starting  and  ending  conditions  is  shown  for 

k+1  ..k+1  >  /l  .  \  .  ,  d  k+1  \  , 


k+I  k+1  \  *  * 

■^N, MS-1/2  ’  UN,ME+l/2  )’  Pn,MS  *  "N,ME+1 

fables  IV-2,  IV-3,  and  IV-4,  respectively. 


,  and 


N,MS  ’  N ,ME+l/2 


Boundary  conditions 


163.  For  a  solid  boundary  at  the  lower  end  of  a  computational  segment 

k+1  k+1 

u  MS  ^  *  0  .  For  a  river  or  other  flow  input,  uM  MC_, /0  =  Q(t)  ,  a 


N, MS- 1/2 


Table  IV- 2 


Table  IV-4 


^ridiagonal  Matrix  Solution  Procedure  (x  - 


"*,MS 


k+l 

“n, ME- 1/2 


Same  as  Step  1  in  Table  IV-3  for  R^,  and  S 


Same  as  Step  2  in  Table  IV- 1  for  m  *  MS+L,  ME-1,  1 
and  as  computed  in  Step  2  Table  IV-2 


k+l 


N  ,ME  ’  ^E^.ME+l/a  +  QME 


k+l 


un ,m+l / 


0  -  -R  , ,  +  S 
2  m  n , m+ 1  m 


k+l 

n,m  PmUn,m+l/2  + 

u’"+l  =  — R  +  S 

N.MS+l/2  ^IS  N.MS+l  MS 


m  =  ME-1,  MS+1 
n  *  N 


specified  flow  velocity.  Alternatively,  the  user  may  specify  the  flow  rate 
and  the  velocity  is  computed  by  dividing  this  by  the  effective  flow  area. 

164.  For  a  tidal  boundary  at  the  lower  end  of  a  computational  segment, 

ys  =  cos  ,  w^  ,  and  ^  are  specified.  Al¬ 

ternatively  the  user  may  specify  a  table  of  values  over  time.  In  addition,  it 
is  possible  to  save  elevations  at  appropriate  locations  on  a  global  grid  and 
later  use  these  elevations  to  drive  a  refined  grid. 

165.  The  treatment  of  the  upper  computational  segment  boundary  is  simi¬ 
lar  to  that  just  presented. 

166.  After  completing  the  above  procedures  for  n  =  2  ,  NMAX  -  1  ,  the 
grid  is  sweeped  in  the  other  direction  (row  wise)  for  each  m  .  In  the 

i,,  -  v  sweep  Equation  set  (1.1.18)  is  considered.  The  computational  stencils 
for  the  continuity  and  momentum  equations  are  as  shown  in  Figures  IV-5, 

and  IV-6,  respectively.  Observe  in  both  Figures  IV-5  and  IV-6,  that  for  a 
given  n  ,  the  variables  at  the  most  advanced  time  level  all  lie  on  a  straight 
line  of  constant  m  .  Since  at  the  most  advanced  time  level,  the  variable  of 


interest  (either 


n+1 72, m 


is  a  function  of  its  nearest  neighbors 


,  the  difference 


along  constant  m  either  /  v  or  n  ,  n  \,  the  difference 

\  n_  1/  ^  ,m  n  ,m  n  >  1  *  m  J 

schemes  are  implicit.  v  J 

167.  The  solution  procedure  for  the  y  -  a7  sweep  is  for  each  m  , 
determine  the  number  of  computational  segments  and  their  beginning  and  end 


locations  in  the  grid.  The  number  of  computational  segments  and  their  loca¬ 
tions  are  a  function  of  the  grid  and  boundary  conditions.  The  computational 
procedures  for  the  y  -  a7  sweep  are  completely  analogous  to  these  presented 
in  Tables  IV- 1  through  IV-4  for  the  x  -  3^  sweep.  Boundary  condition  con¬ 
siderations  also  remain  analogous. 

Programming  considerations 

168.  The  above  algorithms  have  been  programmed  in  the  FORTRAN  language. 
Since  FORTRAN  array  element  variables  cannot  be  referenced  using  fractional 
indices,  the  index  system  shown  in  Figure  IV-7  below  has  been  developed.  Thus 
variables  defined  at  cell  faces  will  take  the  same  index  as  the  center  cell 


variable  to  the  left  or  below,  whichever  is  appropriate.  A  separate  index 
system  is  developed  for  the  stretching  functions  ,  and  _  ,  as  follows: 


XMU(2M-1)  =  (u.) 

1  m 

XMU(2M)  5  ('V„+l/2 


m  and  M  •:  (  1  ,  MMAX 1 


YNU(2N-1)  =  (u,) 

1  n 

n  and  N  •:  (  1  ,  '."MAX) 

YNU(2N)  (^)a+u2 

169.  Thus  the  stretching  function  at  the  center  of  the  cells  has  an  odd 
array  reference,  while  on  the  cell  faces,  its  reference  is  even.  The  follow¬ 
ing  variable  associations  are  employed  for  the  -  x  and  a?  -  v  sweeps: 

-  x  sweep 


(N,M) 

*-*  II 

XU  «-  XMU ( 2*M) 

(N.M+l) 

—  JJ 

XU1  XMU  ( 2*M- 1 ) 

n:+i ,m) 

T  T  2 

YU  ~  YNT ( 2*N- 1 ) 

n:-i,y.) 

*"*  III 

A(M) 

‘  A 

m 

TMP3(M)  =  a^1/2 

TMPl(M) 

am-l/2 

TMP4 (M)  =  a  , 

m+ 1 

yv  p  2  !  M  ) 

am^l/2 

B(M)  =  B 

m+ 1  /  2 

P(M)  ,  Q(M)  ,  R (M)  ,  S (Ml  =  P  ,  0  ,  R  ,  S 


y  sweep 


2 


(N,M) 

«—  JJ 

XU 

XMU ( 2*M- 1 ) 

(N+l.M) 

—  JJ2 

YU1 

YNU  (  2*N-1 ) 

(S-l.M) 

--  JJ1 

YU 

<-*  YNU ( 2*N) 

(N.M+l) 

«->■  LL1 

(N’.M-l) 

LL2 

A(S) 

c 

c 

111 

™P3(N»  s  Vl/2 

TMPl(N) 

s  an-l/2 

TMP4(K)  5  an+1. 

TMP2(N) 

"  an+l/2 

B(N>  E  Vl/2 

P(N),  Q(N),  R(N),  S(N)  E  P  Q  ,  R 

n  n  n  n 


Cell  face  flag  conventions 

170.  In  order  to  control  the  beginning  and  end  of  the  x  -  and 
y  -  sweeps  each  cell  face  is  given  a  two-digit  code.  The  code  is  used 
also  to  determine  boundary  types  and  convective  acceleration  (advection)  ap¬ 
proximations  in  the  vicinity  of  flow  boundaries.  The  code  conventions  are  as 
indicated  in  Table  IV-5  below.  The  u-face  cell  code  is  associated  with  the 

v-velocity  and  is  stored  in  array  ICU:  e.g.  ICU  =  u  .  Similarlv,  the 

n,m  n,m 

v-face  cell  code  is  associated  with  the  v-velocity  and  is  stored  in  array  ICV; 


e.g. 


ICV 


n  ,m 


V 

n  ,m 


Numerical  approximation  of 

the  convective  acceleration  terms 

171.  If  one  considers  Equation  (1.1.14)  for  cij  -  x  sweep  and  Equation 
(1.1.16)  for  the  -  y  sweep,  the  terms  underlined  designated  by  ©,  con¬ 
stitute  the  convective  acceleration  terms.  These  terms  are  approximated  using 


Table  IV- 5 

Face  Ce  1  L  Flag  Convention 


r“-> 


Exposed 

Harrier 


2  Overtopping 
Barrier 


Global 
Grid  Cond  . 
on  LtJM  24 


p?  barrier  no 
'2  ' 


If  n^1)  n,  0 


barrier  no. 

0 


n2 

If  n9  >  9  n-. 


n0  not  used 


(Advec  tion  code) 


Submerged 

Barrier 


n7  barrier  no. 


s  Land  / 

(Water  or  land 
Boundarv  ) 


b  .-.ater 


ADV  =  1 


AOV 


l! 


n->  >  9  n-> 


0  No  advec  tion 

1  nx-direction 

2  ny-d irec  t ion 
*3  nx,  v 

4  nv  ax 

5  nx  ay 

6  ax 

7  av 

8  ax  av 


a  =  anprox. 
n  =  normal 


7  Transmission 
Boundarv 


8  Elevation 


Input 


9_  Flow 
Input 


o  lower  boundary 
1  upper  boundary 

6  Channel  Match  co 

7  Ponding  Condition 

8  K/V 


9  Extrapolated  wave 


n-5  input  no .  n2 


n >  =  0  P  ,  on  bdy 

<_  a 


n2  input  no .  n->  ■  0 

n0  =  Q  =  o  on  b 
2  Bn 


1)4 


centered  dit f erencing.  However,  in  the  vicinity  of  external  boundaries 
("land 'water  interface)  or  near  internal  barriers,  either  the  variables  are  not 
defined  at  appropriate  locations  on  the  computational  stencils  Fig*  IV-3-lV-n 
or  it  may  be  inappropriate  to  form  the  differences  across  the  flow  obstruc¬ 
tions.  In  the  hydrodynamic  program,  these  terms  are  not  approximated  in  such 
cases;  that  is,  the  motion  equations  are  linearized.  In  addition,  the  lateral 
diffusion  terms  (the  underlined  terms  designated  by  G)  )  are  also  ignored, 
when  the  convective  acceleration  terms  are  not  approximated.  This  is  accom¬ 
plished  by  specifying  ADV  =  I  on  input.  Subroutine  ADVBAR  is  used  to  set 
the  n0  digit  for  wet  cells  surrounding  barriers. 

172.  One-sided  difference  approximations  are  made  for  both  convective 
acceleration  and  lateral  diffusion  terms  if  the  user  specifies  ADV  =  2  . 

These  options  are  still  in  the  test  and  research  phase  and  are  not  intended  at 
this  time  for  general  use. 

Computational  line  development 

173.  Let  us  first  consider  the  form  of  the  recursion  equations  used  in 
the  program  for  the  cx^  -  x  sweep.  Using  the  program  variable  notation: 


SEP(N.M)  =  -P (M)  *  UP (N ,M)  +  Q(M) 


UP(N.M-l)  =  -R(M-L)  *  SEP*(N',M)  +  S(M-l) 


M  =  ME,  MS,  -1 


(1.7.1) 


174.  Consider  the  general  form  of  a  computational  segment  as  shown  in 

Figure  IV-8.  Let  us  now  consider  the  case  of  a  solid  boundary  at  the  lower 

end  of  the  computational  segment.  Since  the  ICU(N,M)  flag  is  searched  until 

ICU CN,M) / 10  equals  6  ,  MS  =  M  ,  such  that  ICU(N,M)/10  is  6  .  In  this 

case,  UP(N,MSS)  =  0  .  In  order  to  achieve  this,  R(MSS)  =  S(MSS)  =  0  .  This 

k+1 

corresponds  to  the  case  in  which  u  ,  =  0  .  Consider  the  case,  when  the 

r  n, ms- 1/2 

user  specifies  a  velocity,  u*  ,  at  u-cell  face  (N,M).  In  this  instance, 

MS  =  MS+L  ,  and  R(MSS)  =  0  ,  S(MSS)  =  u*  where  MSS  =  MS-1  .  This  cor- 

k+I 

responds  to  the  case  for  which  u  , ,  is  specified.  Alternatively,  con- 
sider  the  case  when  the  user  specifies  an  elevation,  n*  ,  at  the  center  of 
cell  IN, MS)  in  this  case  MS  »  MS+1  and  R(MS-l)  and  S(MS-l)  are  computed 
directly  in  terms  of  "*  .  This  corresponds  to  the  case  in  which  "*  is 


s  (MSS)  Q  (MS)  S  IMS)  Q  (ME)  S  (ME)  Q  (MEE) 

MS  -  start  of  the  computational  segment 
ME  =  END  OF  THE  COMPUTATIONAL  SEGMENT 

Figure  IV-8.  General  form  of  a  computational  segment 
for  the  x  -  sweep 

175.  Consider  now  the  upper  end  of  the  computational  boundary.  In  the 

case  of  a  solid  boundary,  ME  =  M  ,  such  that  ICU(N,M)/10  equals  5  .  In 

this  case,  l'P(N,ME)  =  0  .  Notice,  however,  in  the  above  equation  pair, 

t’P(N,M-l)  and  SEP(N,M)  ,  correspond  to  the  center  and  left  side  of  cell 

(N,M).  Therefore,  we  may  compute  SEP(N,ME)  directly  in  terms  of  the  zero 

boundary  velocity.  The  specification  of  a  flow  velocity  is  treated  in  exactly 

the  same  manner,  since  it  is  specified  at  the  same  location  as  a  solid  bound- 

k+L 

ary.  Both  these  cases  correspond  to  specifying  u^  rae+i/?  •  Consider  the 

specification  of  an  elevation,  n*  ,  at  cell  (N,ME).  In  this  again, 

ME  =  ME-1  ,  and  it  is  necessary  to  first  calculate,  UP(N,ME)  in  terms  of 

'*  .  This  corresponds  to  specifying  n*  .  Observe  analogous  considera- 

n,me 

tions  hold  for  the  -  y  sweep. 

Sub-grlj  barriers 

176.  Butler  [34]  has  represented  flow  restrictions  on  cell  faces  through 
a  thin  wall  barrier  concept.  The  width  of  the  barrier  essentially  is  assumed 
zero  and  does  not  enter  the  computations.  Barriers  in  the  Lake  Okeechobee 
Study  are  of  the  following  two  types: 

a.  Exposed  barrier. 

b.  Dynamic  barrier. 


136 


Exposed  barriers  are  considered  to  fall  within  two  categories.  Category  one 
barriers  r.av  be  specified  to  breech  in  a  manner  presented  subsequently.  Cati 
gory  two  barriers  may  be  overtopped  without  breeching.  Flow  information  for 
both  categories  is  written  to  tape  for  subsequent  input  to  a  link-node  drain' 
age  model.  Dynamic  barriers  may  alternate  between  submerged  and  exposed  con' 
ditions  through  the  simulation.  Within  the  Lake  Okeechobee  Study,  these  bar' 
riers  are  used  to  represent  the  tree  islands  and  are  assumed  initially  to  be 
exposed . 

177.  The  model  user  specifies  barrier  face  condition  through  input  as 
follows : 

ITYP  e  Type  of  barrier  face  (1  exposed,  2  dynamic) 

IDIR  =  Direction  of  restriction  (1  flow  restricted  on  u-face,  2  flow 
restricted  on  v-face 

INTDX  E  Index  number  used  to  specify  the  admittance  coefficient  for 
submerged  conditions 

II  =  Grid  index  of  barrier  line 

12  I 

)  E  Barrier  extends  from  cell  12  to  13 

13  I 

14  e  Barrier  configuration  index 

0  link-node  overtopping  barrier 

15  E 

NCR  breechable  link-node  barrier  with  NCR  cross-sections 

16  E  Barrier  face  condition  and  transmission  code 
DUM1  =  Barrier  toe  elevation  (rasl) 

DUM2  H  Barrier  lakeside  side  slope 

DUM3  E  Nearshore  beach  slope 

DUT14  E  Barrier  crest  elevation  (msl) 

The  above  information  is  saved  in  the  following  arrays  for  each  barrier  face 
where  KST  is  a  barrier  face  counter: 

I3ARR (KBT ,  1)  <-  IDIR 
IBARR(KBT , 2)  -  ITYP 
I  BARR  ( KBT  ,  3)  <-  INDX 
IBARR  (KBT  ,4)  <-  N 
IBARR(KBT , 5)  -  M 

IBARR(KBT ,6)  -  16  (If  16  =  0,  IBARR(KBT ,6)  =  17) 

IBARR(KBT,  7)  «-  15 

KBARO (KBT)  -  14 


BARR (KBT , 1 )  -  DUM1 
3ARR(KBT,2)  -  DUM2 

BARR  i' KBT  ,3)  -  CVMGZ  (DUM2  ,  DUM3  ,  DUM3) 


Ol 


All  type  1  barriers  are  either  overtopped  or  breechable  barriers 

KBAROO(KGG)  <-  KBT  Link-node  model  overtopping  barrier  index 
-  .  KBAR(KFTB)  +-  KBT  Link-node  model  breechable  barrier  index 
For  each  cross-section  J=L  ,  15  ,  the  user  specifies: 

(ELEVCR( J ,KFTB) ,  J=l,  15)  the  final  cross-section  elevation  after 
breeching 

178.  A  cell  face  flag  is  computed  as  follows: 


IDIR  =  1,  u  face  barrier 
ICIL.  w  =  10  *  ITYP  +  Ih’DX 


IDIR  =  2,  v  face  barrier 

ICV,  w  «  10  *  ITYP  +  IN'DX 
N  ,M 


179.  For  the  case  in  which  advection  (the  convective  terms  in  the 
motion  equation)  is  considered,  one  may  employ  subroutine  ADVBAR  to  linearize 
the  motion  equation  approximation  in  the  vicinity  of  the  barriers.  We  present 
in  turn  the  computational  procedures  for  each  barrier  face. 

Exposed  Barrier  Face 

180.  This  barrier  face  may  never  experience  any  type  of  flow.  The  face 

flag  (ICU  or  ICV  )  equals  In  ,  where  r  <  (0,9)  .  This  face  flag  is 
n  f  rn  n  >  in 

divided  by  10  to  yield  a  face  type  of  1.  The  face  type  is  examined  in  the  x 
and  y  sweeps  of  the  grid  within  the  ADI  procedure  in  order  to  establish  compu¬ 
tational  segments.  For  a  face  type  of  1  the  computational  segment  is  either 

k+ 1  k+ 1 

initiated  or  terminated  with  u  *  0  or  v  =  0  .  The  velocity  at  bar- 

n  f  in  n  ,m 

rier  face  is  thus  set  to  zero.  This  condition  is  maintained  throughout  the 
simulation.  The  barrier  thus  allows  no  flow  (velocity  zero)  through  itself. 

LSI.  The  barriers  may  also  be  breeched.  The  model  user  specifies  an 
initial  breech  time  In  hours  relative  to  the  start  of  the  simulation  and  a 

miration  ■  f  breeching  also  in  hours.  For  each  cross-section  in  each  hreech- 
ib  lu  barrier  'tie  computes  a  breech  rate  per  time  step.  This  rite  is  assumed 
•.  he  ur.i:  rm  over  the  duration  of  breeching  activity. 

132.  In  this  manner  it  is  possible  to  simulate  breech  character ist ics 
.arming  -ver  manv  grid  cell  faces  or  on  1 v  over  a  nart  of  a  single  grid  cell 


tace.  The  types  of  failure  cross-sections  that  can  be  simulated  are  shown  in 


Figure  IV-9.  The  following  notation  is  used  in  the  figure,  h*.  ,  corresponds 
to  the  levee  height  at  time  k  for  levee  i  for  cross-sect  ion  j  . 

133.  Each  levee  section  associated  with  a  grid  cell  face  will  be 
associated  with  a  node  in  the  link-node  model.  The  flow  rate  (I.^'T)  and  water 
surface  elevation  associated  with  either  overtopping  or  breeching  will  be 
written  to  disc  for  each  levee  section  along  a  grid  cell  face  at  each  time 
step  in  the  hydrodynamic  simulation  of  the  lake.  The  flow  rate  computations 
are  as  follows: 


Multi-cell  face  breech: 


q 


k 

i 


(l//T) 


(a) 


k  k£  .  3  ,  . 
v i  =  q1  (~  /T) 


(b) 


with  h.  .  =  TB 


(1.8. L) 


BREECH 


LEVEE  CROSS-SECTIONS 


C0R  ALL 


Figure  IV-9 .  Cell  breech  mechanics 


140 


m 


Single  cell  face  breech: 


with  h~  =  ZBi 
where 


<i  ■  -  < 


vk  =  qk'S  (L3/T) 
ill 


(L  /T) 


(1.8.2) 


=  flow  rate  per  unit  width  of  breech  at  time  step  k  for  cell  i 

"  =  water  surface  elevation  at  time  step  k  for  cell  i 

ZB^  =  levee  elevation  for  cell  i  (msl) 

ZB^  =  levee  elevation  for  cell  i  cross  section  increment  j  (msl) 

nc .  =  number  of  increments  in  cross  section  for  cell  i 
1 

k 

h  =*  breech  height  for  cell  i  at  time  step  k 

h  .  =  breech  height  for  cell  i  at  time  step  k  cross  section 

^  increment  j 


2-S^  =  cell  face  width  along  the  breech  at  cell  i 

v^  =  flow  rate  caused  by  breeching  during  time  step  k  at  cell  i 

The  link  node  model  will  access  the  above  disc  file.  The  inflow  at  an 
arbitrary  node  I  from  Lake  Okeechobee  will  be  developed  as  the  sum  of  all 
flows  at  VIFM  cell  faces  assigned  to  node  I.  Thus  assuming  equal  time 
steos  '300  sec)  in  the  WIFM  and  link-node  models. 


i 

£Vi 


k  <  hk 
I  i 


n  k  ,  k 

0  uI-hi 
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where 


Qj.  =  flow  rate  from  Lake  Okeechobee  into  node  I  at  tine  step  k 
'K 

r.  =  factor  to  either  allow  or  disallow  flow  from  cell  face  i  at  time 
step  k  into  node  I 

k 

=  flow  rate  at  cell  face  i  at  time  step  k 

h  -  elevation  at  ceil  face  i  at  time  step  k 

k 

H  =  N'ode  I  elevation  at  time  step  k 

As  may  be  observed  from  the  above  equation,  flow  from  the  lake  at  cell  face  i 

at  time  step  k  is  not  allowed  if  the  node  elevation  is  greater  than  the  cell 

elevation  in  WIFM.  The  coupling  mechanism  is  essentially  one-way  with  water 
passing  from  the  lake  to  the  link-node  drainage  model.  Water  is  not  allowed 
to  pass  from  the  drainage  model  into  the  lake.  This  coupling  provides  a  con¬ 
servative  estimate  of  flood  levels  outside  the  lake  appropriate  for  evacuation 
and  water  conservation  area  planning. 


Dynamic  Barrier  Face 


184.  This  type  of  barrier  face  changes  dynamically  as  water  levels  rise 
and  fall  throughout  the  course  of  the  simulation.  A  counting  mechanism  is 
employed  such  that  a  given  face  type  must  hold  for  a  user  specified  number  of 
time  steps.  This  procedure  is  similar  to  those  used  in  the  flood/drv  scheme 
and  eliminates  instabilies  in  the  solution  introduced  by  discrete  changes  in 
barrier  face  type. 

185.  The  computational  procedures  are  performed  in  Subroutine  FLOODY 
prior  to  considering  the  flood/dry  cells.  The  order  in  which  the  barriers  are 
processed  is  reversed  each  time  step  to  remove  computational  bias.  All  bar¬ 
riers  are  indexed  within  the  DO  650  loop.  For  ITYP  =  1  (exposed  barriers  in 
a  breechable  state)  control  is  transferred  to  the  end  of  the  loop  and  the  next 
barrier  is  considered.  These  barriers  are  handled  prior  to  entering  the 

DO  650  loop  as  discussed  previously.  Dynamic  barriers  are  in  one  of  two  pos¬ 
sible  states  (ITYP  =  2  exposed  condition  or  ITYP  =  4  submerged). 

186.  The  following  elevations  are  defined  in  order  to  develop  the 
methodology  for  handling  these  barriers: 
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■1V1 


k+1 

T3  =  n  ,  +  Aa  , .  (c) 

n+l,m  n+l,m 


T4  "  hB  +  S 


(d) 


where 


k+l 

n 

n  #m 


water  surface  elevation  at  location  (n,m)  at  time  k+1 
respect  to  model  datum  (12.5  msl) 


Aa 

n  ,m 


cumulative  water  surface  adjustments  at  location  (n,m) 
k+1 


dt 


h 


B 


'b 

Next  define 


=  model  datum  adjustment  to  barrier  datum  (12.5) 

=  carrier  elevation  (msl) 

=  minimum  amount  of  water  for  water  surface  adjustments 
the  following  variables: 


Case  I:  IDIR  =  1  (Barrier  to  the  x  -  velocity  component) 


T5  -  T2  +  dt 


T6  «-  ICU 


IG  - 


II 

( n  ,  in+ 1 ) 


DS1 


Aa,  (u  ) 

1  1  m 
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with 

at  time 
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Case  II:  IDIR 


(Barrier  to  the 


v  - 


velocity  component) 


T5  -  T3  +  dt 


T6  -  ICV 


IG 


II 
(n+1 ,n) 


DSI 


i-s-  (u0) 

2  2  n 


DS2 


Aa, (u0) 
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DL  «-  ia.  (u, ) 

I  1  m 

Let  us  now  consider  the  submerged  condition  (ITYP  =  4)  .  IBART^.  is  incre¬ 
mented  and  compared  to  KS3  ,  which  is  specified  through  user  input.  If 

IBART^.  =  KS3  ,  then  a  check  on  T1  and  T5  is  performed  else  no  further  pro¬ 

cessing  is  performed  on  this  barrier.  If  both  TI,  T5  >  T4  ,  the  barrier  is 
still  submerged  and  the  next  barrier  is  considered. 

187.  If  the  above  condition  does  not  hold,  the  barrier  type  is  set  to  2 

corresponding  to  an  overtopped  condition.  The  face  change  counter  IBART  is 

«et  to  one.  The  cell  face  flag  is  also  appropriately  modified;  e.g.,  20  is 

subtracted  from  ICU  or  ICV  .  The  next  barrier  is  then  considered. 

rt  y  in  n,m 

188.  Consider  next  the  overtopped  condition  (ITYP  =  2).  IBART^  is 

incremented  and  compared  to  KS4  ,  which  is  specified  through  user  input.  If 

I3ARTt  =  KS4  ,  then  a  check  on  TL  and  T5  is  performed  else  no  further  pro¬ 

cessing  is  performed  on  this  barrier.  If  both  Tl  and  T5  are  greater  than 
74  ,  the  barrier  condition  is  changed  to  submerged  and  the  IBART^.  counter 
reset  to  I  .  For  TI  or  T5  <  T4  ,  the  following  quantities  are  computed: 


h  =  MX  (TI,  T5) 
sax 


h  <  h„  ,  consider  the  next  barrier 
sax  3 


h  .  =  min  (Tl ,  T5) 

run 


If  h,  <  h  ,  then  set  h  =  0.5  *  (T1  +  T2) 
b  min  max 
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dbc  =  h 


max 


3/2 

g  =  CO  dbc  ,  which  corresponds  to  the  flow  per  unit 
L  NUA 

barrier  length 


Ar  =  gx/DSl,  ^here  t  is  the  time  step  length 


=  -gt/DS2 


The  elevation  increments  are  adjusted  in  a  cumulative  fashion  as  follows: 


Aa  =  Aa  +  A",  (a) 

n ,  m  n  ,  m  1 


IG 


*  Aa 


IG 
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The  most  advanced  time  level  velocities  UPTT  or  VP  are  set  to  zero.  The 

depth  adjustments  An-^  and  A",,  are  limited  such  that  the  elevation  adjustment 

the  cell  with  the  largest  water  surface  elevation,  will  not  cause  that  elevation 

to  be  less  than  h,  +  c,  . 

b  b 

flooding  and  drying  methodology 

L89.  Leendertse  [35]  was  perhaps  the  first  investigator  to  develop  a 
flooding/drying  capability  in  the  simulation  of  a  tidal  flat  problem.  The 
initial  approach  was  to  make  discrete  changes  (on  the  faces  of  grid  cells'* 
between  land  and  water  directly  on  the  finite  difference  grid.  In  so  doing, 
one  must  be  aware  that  the  discrete  changes  in  the  computational  boundaries 
may  introduce  noise  into  the  finite  difference  computation.  In  general,  when 
new  cells  enter  the  computation,  their  initial  water  depth  is  relatively  small 
resulting  in  a  large  value  of  bottom  friction,  which  in  turn  will  dampen  out 
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oscillations  < noise)  introduced  bv  the  new  boundary  change.  This  frictional 
dampening  mechanism  is  in  agreement  with  the  physical  process  and  does  natu¬ 
rally  stabilize  the  computations.  However,  in  order  to  allow  this  mechanism 
time  to  work,  it  is  usually  necessary  not  to  change  the  boundarv  too  often; 
i.e.,  not  every  time  step.  This  approach  is  also  in  keeping  with  the  general 
philosophy  of  applying  finite  difference  techniques  to  relatively  slowlv 
changing  processes. 

190.  In  developing  a  technique  for  treating  a  moving  boundarv  problem, 
the  above  concepts  are  assumed  to  form  the  computational  foundation.  In  addi 
tion,  the  consideration  of  economy  must  also  be  introduced,  since  the 
flooding/ drying  techniques  must  be  incorporated  in  addition  to  the  finite 
difference  computations  over  the  interior.  In  order  to  formulate  the  most 
appropriate  approach  for  Lake  Okeechobee,  a  review  of  techniques  is  first  pre 
ser.ted.  The  details  of  the  Lake  Okeechobee  problem  are  next  presented  fol¬ 
lowed  by  the  development  of  the  most  appropriate  moving  boundary  techniques. 


191.  The  search  for  the  next  boundary  is  made  at  intervals  larger  than 
the  time  step.  This  enables  the  computational  noise  generated  by  the  boundar 
change  to  die  off  before  it  is  time  to  reexamine  the  boundary  location.  Thus 
the  selection  of  the  next  boundary  is  not  influenced  by  the  perturbations 
caused  by  the  last  boundary  change. 

.Q2.  If  in  a  particular  field  (cell)  a  cross  section  decreases  to  less 
than  a  preset  value  (greater  than  zero),  then  that  field  (cell)  is  removed 
from  the  computation.  During  rising  water  levels,  grid  points  are  added  to 
the  computation  if  the  average  of  the  adjacent  fields  which  are  underwater 
results  for  ail  four  sides  of  the  field  (cell)  in  cross  sections  which  are 
larger  than  the  above-mentioned  preset  value.  If  this  is  the  case,  the  water 

level  in  the  new  field  (cell)  is  set  to  the  average  value  of  the  adjacent 
fields  fcells)  in  the  computation. 

I'5?.  The  principle  of  these  tvpe  of  boundary  changes  while  very  simple 
to  write  is  very  complex  to  implement  computationally. 


194.  In  order  to  aid  the  description  of  the  revised  land-water  bound¬ 
aries,  a  representation  portion  of  the  space  staggered  grid  employed  is  shown 


in  Figure  IV- 10. 

195.  The  volume  element  associated  with  each  grid  cell  j,k  at  time 
level  n  is  defined  as  follows: 


n 


v°v,k  ■  r- . 


<c  +  4  (hj-l/2,k-l/2  +  hj-l/2,k+i/2  +  hj+l/2,k-I/2 


+  hj+l/2,k+l/2) 'Ax 


(1.9.1) 


The  cross  section  depth  between  grid  cell  j ,k  and  j  +  l,k  at  time  level  n 
is  given  by  the  following  relation: 
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+  n 
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j+l/2,k-l/2^ 


(1.9.2) 


The  three  cross  section  depths  between  grid  cell  j,k  and  the  other  three 
surrounding  grid  cells  are  defined  similarly.  The  Chezv  coefficient,  C  , 
which  is  computed  at  each  water-level  location,  is  used  to  designate  whether  a 
grid  cell  is  wet  or  dry.  A  value  of  E  =  0  designates  land,  while  a  positive 
value  signifies  water. 

196.  There  are  two  checks  made  to  determine  whether  a  grid  cell  that  is 
currently  wet  should  become  dry.  If  after  each  calculation  of  the  new  water 
level  a  negative  volume  is  obtained  for  the  grid  cell,  then  that  cell  is  re¬ 
moved  from  the  computation  (made-  dry) .  This  is  performed  by  setting  the  Chezv 
coefficient  to  zero  as  well  as  the  velocity  components  on  each  of  the  four 
cell  faces.  In  order  to  include  the  influence  of  the  new  land  point  on  the 
adjacent  water  levels  and  velocities,  the  previous  row  or  column  (note  an  ADI 
scheme  is  being  employed)  of  recursion  formulas,  velocities,  and  water  levels 
is  recalculated  with  the  point  taken  as  dry.  In  addition,  the  calculation  of 
all  velocities  and  water  levels  along  present  row  or  column  (row  or  column  in 
which  the  grid  cell  became  dry)  are  repeated  before  the  computation  proceeds 
to  the  next  row  or  column.  When  a  grid  cell  becomes  dry,  it  is  assumed  that  a 
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Figure  IV-10.  Representative  Space  Staggered 
Grid  System 
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thin  layer  of  water  will  remain  on  it.  The  depth  of  this  water  is  set  to  the 
level  calculated  in  the  previous  time  step,  since  this  calculation  did  not 
vield  a  negative  volume  element.  The  second  check  employes  the  cross  section 
depth  concept.  If  any  one  of  the  four  cross  section  depths  associated  with  a 
grid  cell  become  negative,  the  grid  cell  is  taken  out  of  the  computations; 
i.e.,  the  C  value  at  this  grid  cell  is  set  to  zero  along  with  all  cell  face 
velocities. 

197.  Since  a  negative  volume  or  negative  cross  section  might  be  calcu¬ 
lated  at  any  point  in  time  in  the  computation,  the  above  two  checks  are  made 
at  each  time  step.  Therefore  discrete  changes  in  the  system  boundaries  can 
occur  at  each  time  step.  In  order  to  prevent  the  occurrence  of  a  large  number 
of  changes  in  consecutive  time  steps,  the  following  boundary  change  patterned 
after  the  original  procedure  is  used. 

198.  If  any  of  the  four  cross  section  depths  associated  with  a  grid 
cell  become  less  than  a  positive  preset  value,  then  that  grid  point  is  taken 
out  of  the  computation.  Since  a  positive  preset  value  is  used,  this  condition 
is  more  restrictive  than  the  above  two  checks.  The  water  level  is  set  at  its 
current  value  and  is  maintained  at  this  level  until  the  point  floods.  By  per¬ 
forming  this  more  restrictive  procedure  at  intervals  larger  than  the  time 
step,  it  is  found  that  a  large  proportion  of  the  boundary  changes  are  made, 
while  only  a  few  changes  occur  because  of  negative  volumes  or  cross  section 
depths.  Thus  the  effect  of  local  discontinuities  generated  by  discrete  bound¬ 
ary  changes  is  reduced. 

199.  The  flooding  procedure  also  was  revised  from  the  original  method. 
For  a  dry  grid  cell  each  of  the  four  surrounding  grid  cells  are  checked  to  see 
if  they  are  wet.  If  one  or  more  of  the  neighboring  cells  are  wet,  then  the 
water  levels  of  these  grid  cells  are  averaged.  If  this  average  water  level  is 
larger  than  the  level  retained  on  the  dry  grid  cell,  then  this  grid  cell  may 
be  allowed  to  flood.  A  check  is  made  on  the  cross  section  depth  between  the 
dry  grid  cell  and  each  of  the  surrounding  wet  grid  cells.  If  any  of  these 
cross  sections  are  negative,  the  point  remains  dry.  If  all  are  positive,  the 
grid  cell  Is  added  back  to  the  computational  field.  The  water  level  is  set  at 
the  value  which  was  allowed  to  remain  when  the  cell  became  dry.  Note  it 
appears  that  a  cell  must  dry  before  it  can  flood;  i.e.,  the  tidal  flat 
simulation  is  initiated  at  high  tide.  The  check  for  flooding  is  also  made  at 


intervals  larger  than  the  time  step  in  order  to  reduce  the  discrete  boundary 
change  effects  introduced  by  new  flood  cells. 


Original  Butler  Technique  [34] 

100.  In  the  development  of  the  Waterways  Implicit  Flooding  Model 
(WIFM) ,  an  alternating  direction  implicit  solution  technique  has  been  em¬ 
ployed.  A  cell  face  flag  system  has  been  developed  in  order  to  determine  in 
each  sweep  where  the  computational  line  should  begin  and  end.  This  cell  face 
flag  system  is  modified  in  the  flood/dry  procedure  to  designate  the  moving 
boundary.  Once  a  cell  face  has  been  opened  it  may  not  be  closed  for  a  user 
specified  number  of  time  steps  (greater  than  one).  Similarly,  if  a  face  has 
been  closed  it  may  not  be  opened  for  a  user  specified  number  of  time  steps 
(greater  than  one).  These  procedures  are  similar  to  those  developed  by 
Leendertse  [35,36]  and  tend  to  reduce  the  noise  generated  by  discrete  changes  in 
the  boundary  affected  by  opening  and  closing  cell  faces. 

20 1 .  Butler  [34]  has  associated  four  cell  face  states.  A  face  may  be 
partially  open,  open,  partially  closed  or  closed.  In  the  partially  open 
state,  a  weir  formula  is  used  to  effect  a  volume  transfer  between  the  adjoin¬ 
ing  cells.  No  cell  face  velocity  is  computed.  In  the  open  state,  the  appro¬ 
priate  momentum  equation  is  used  to  determine  the  cell  face  velocity.  In  the 
partially  closed  state,  a  volume  transfer  concept  based  on  a  fraction  of  cell 
water  depth  is  used  to  drain  water  between  adjoining  cells.  No  cell  face 
velocity  is  computed.  In  the  closed  state,  the  cell  face  velocity  is  zero. 
L'nder  this  approach  a  cell  face  velocity  is  associated  only  when  a  cell  face 
is  open  or  closed. 

202.  As  noted  above,  a  cell  face  concept  is  used  in  implementing  the 
computations.  For  a  given  cell  (i,j),  a  u  face  (i,j+l/2)  and  a  v  face 
(i+I/2,j)  are  designated  in  the  code  with  the  same  (I ,J)  index.  Thus  one  may 
associate  each  cell  (I,J)  with  two  faces.  Butler  [34]  sets  up  a  list  of  cells 
which  may  either  flood  or  dry.  In  effect,  these  cells  are  used  to  identify 
the  faces  which  may  change  state  during  the  simulation. 

203.  The  model  user  specifies  the  following  variables  to  control  the 
flood/dry  operation: 
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EPSD  =  minimum  amount  of  water  allowed  on  a  cell  in  order  for 
the  cell  to  still  be  considered  wet 

(Note  this  corresponds  to  Leendertse's  [35]  small 
positive  preset  value) 

XLAND  =  an  elevation  such  that  if  the  bottom  elevation  of  the 

cell  is  greater  than  XLAN'D  ,  the  cell  will  never  flood 

XSCOUR  h  an  elevation  such  that  if  the  bottom  elevation  of  the 
cell  is  smaller  than  XSCOUR  ,  the  cell  will  never  dry 

CD  ,  i=l,20  =  weir  coefficients  for  initial  transfer  of  water  onto  a 
dry  cell  which  is  to  become  wet 

CAYD .  ,  i=L,20  =  recession  coefficients  for  draining  flooded  cells  which 
1  will  become  dry 

KS1  =  minimum  number  of  time  steps  to  hold  a  cell  face  open 
which  had  previously  been  closed 

KS2  z  minimum  number  of  time  steps  to  hold  a  cell  face  closed 
which  had  previously  been  open 

204.  The  implementation  of  the  flood/dry  procedures  is  effected  in  Sub¬ 
routine  FL00DY,  which  is  called  for  each  time  step  in  the  simulation.  The 
general  structure  of  this  subroutine  is  shown  in  Table  IV-6.  Computations  are 
first  performed  for  the  u  cell  faces  and  then  for  the  v  cell  faces.  The 
weir  computations  are  order  independent,  since  the  elevations  used  do  not  ac¬ 
count  for  cumulative  modifications.  The  drain  computations,  however,  are 
order  dependent. 

205.  In  analogy  to  the  Leendertse  [35,36]  procedures,  the  cell  Chezy  C  is 
set  to  zero  for  a  flood  cell,  which  has  become  dry.  Alternatively,  if  a  cell 
face  is  opened  and  constitutes  one  of  the  four  sides  of  a  given  cell,  then 
that  cell  is  considered  wet  and  assigned  a  positive  Chezy  C. 

206.  In  general,  the  history  of  a  cell  face  and  the  hydrodynamic  compu¬ 
tations  are  impossible  to  specify  a  priori.  They  are  an  extremely  complex 
function  of  topography  and  storm  characteristics.  Instances  may  occur  in 
which  it  becomes  difficult  to  distinguish  between  the  partially  open  and  par¬ 
tially  closed  cell  face  states.  In  each  case  u  =  0  ,  ICU  =  50  ;  how- 

J  n  ,m  n ,m 

ever,  one  must  determine  whether  to  use  the  drainage  concept  or  the  weir 
formulation.  Within  the  drainage  concept,  water  is  transferred  downhill  as 
determined  by  topography  and  not  water  surface  elevation.  The  user  specifies 
recession  coefficients  as  previously  mentioned,  which  determine  the  fraction 


Table  IV-6 

Subroutine  FLOODY  General  Structure 


Flood/Dry  Cell  Face  Logic 

A.  Examine  u  Cell  Face  State  (try  to  open) 

(i)  Partially  open — weir  equation 

(ii)  Partially  closed — drain 

(iii)  Completely  open — momentum  equation 

B.  Examine  v  Cell  Face  State  (try  to  open) 

(i)  Partially  open — weir  equation 

(ii)  Partially  closed — drain 

(iii)  Completely  open — momentum  equation 

C.  Examine  u  Cell  face  State  (try  to  close) 

(i)  Partially  closed — drain 
(ii)  Close — set  face  velocity  to  zero 

D.  Examine  v  Cell  Face  State  (try  to  close) 

(i)  Partially  closed — drain 
(ii)  Closed — set  face  velocity  to  zero 

Perform  Flood/Dry  Cell  Chezv  C  Check. 

Update  Surface  Elevations 


of  the  "uphill"  water  depth  to  transfer.  In  general,  the  characterization  of 
rlooding  and  drving  may  be  controlled  by  the  user  through  the  specification  of 
F31  ,  K S 2  ,  FPSD  ,  and  the  CAYD  array.  The  user  also  has  some  control 
over  the  drainage  through  the  specification  of  the  bottom  elevation  of  each 
cell  in  the  flood  plain. 

Revised  Technique  by  Butler  and  Prater  [37] 

20’.  In  conjunction  with  a  surge  study  for  the  southern  shore  of  Long 
Island,  the  flooding  and  drying  methods  have  been  substantially  revised.  The 
major  modifications  are  listed  in  Table  IV-7  below.  Elevation  variables 
instrumented  in  the  computations  have  been  revised  as  shown  in  Figure  IV— 11. 

Table  IV-7 

WIFM  Flood/Dry  Modifications 


The  face  states  are  redefined  as  follows: 

IPS,  JPS  =  0  =  Cell  vulnerable  to  a  change  in  state 

IPS,  JPS  =  1  =  Partially  open  (weir  equation) 

IPS,  JPS  =  2  =  Partially  closed  (drainage  relations) 

The  weir  relations  are  now  used  to  assign  a  cell  face  velocity.  This 

velocity  is  limited  by  assuming  (3u/9t  *  -  g3h/9x  or  3v/3t  =  -  g3h/5y) 

weir  relations  are  now  cumulative. 

The  drainage  relations  are  modified  to  drain  based  upon  water  surface  as 
opposed  to  bottom  cell  elevations.  A  drainage  equalization  check  is 
also  made. 

All  water  surface  elevation  adjustments  are  accumulated  (cumulative)  as 
they  are  performed.  Therefore  the  cell  face  computations  are  no  longer 
order  independent.  Thus  in  alternate  time  steps  u  faces  are  consid¬ 
ered  before  v  faces.  In  addition  the  order  of  the  potential  flood/drv 
cells  is  reversed  sequentially. 
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PART  a:  FACE  OPENING  LOGIC 


PART  b:  FACE  DRAINING  LOGIC 


Figure  IV— 11.  Revised  elevation  concepts  (Continued) 


PART  c:  FACE  CLOSING  LOGIC 


Figure  IV- L I  (Concluded) 


-OH.  .he  general  nature  of  the  computations  is  similar  to  the  original 
procedures.  These  procedures  are  subject  to  user  control  through  the  specif: 
cation  of  ESI  ,  KS2  ,  CD  ,  CAYD  ,  and  EPSD  . 

209.  Butler  [34,37]  is  the  first  investigator  to  our  knowledge  to  have 
employed  a  drainage  concept.  The  CAYD  array  controls  the  fraction  of  the 
water  depth  to  be  transferred  in  one  time  step  and  is  somewhat  empirical  in 
nature.  Several  investigators  [38,39]  have  used  the  weir  equation  to  initially 
transfer  water  onto  a  dry  cell.  This  mechanism  is  used  to  establish  suffi¬ 
cient  water  on  a  dry  cell  in  a  mass  conservative  manner,  such  that  the  momen¬ 
tum  equation  may  be  solved. 

210.  Both  the  original  and  revised  flood/dry  procedures  are  extremely 
complex.  It  is  indeed  difficult  to  account  for  all  flooding  and  drying 
poss ib ilit ies . 

Lake  Okeechobee  Study  Procedures 

211.  Since  the  computations  are  identical  for  the  u  and  v  cell 
faces,  the  detailed  description  will  be  presented  only  for  the  u  cell  face. 
Within  the  DO  45  Loop  in  Subroutine  DIGEST,  the  water  surface  elevations  for 
all  three  time'  levels  are  initialized  to  the  cell  bottom  elevation.  Within 
the  DO  *50  Loop  in  Subroutine  DIGEST  a  flood  elevation  range  is  established 
for  each  cell  face  of  each  potential  flood/dry  cell  as  follows: 


DFLU  =  ZFU  =  max(HTT,  HTT.K_,AV)  +  EPSD  +  0.2 
I  II  II+MMAX 


DFLV  =  ZFV  =  max(HI];,  HII+J)  +  EPSD  +  0.2 


DFLl'F  =  DFLU  -  EPSD  =  tnax(H  ,  HII+J>1AX)  +0.2 


DELVE  -  DFLV  -  EPSD  =  raax(H  ,  H1J+1)  +0.2 


where 


Upper  u  face  elevation  level  for  flooding/drying 
Upper  v  face  elevation  level  for  f looding/drving 


DFLL’E  r  Lower  u  face  elevation  level  for  flooding/drying 
DFLVE  2  Lower  v  face  elevation  level  for  f looding/drving 


Note  an  offset  of  0.2  feet  has  been  used  with  respect  to  the  cell  bottom 
elevation . 

212.  The  cell  face  states  are  saved  in  variable  IPS  for  the  u  face 
and  JPS  for  the  v  face.  These  variables  are  initialized  to  zero  and  are 
updated  at  the  end  of  each  pass  through  Subroutine  FLOODY,  In  addition,  a 
cell  face  counter  (IQU  for  the  u  face  and  IQV  for  the  v  face)  is  ini¬ 
tialized  to  zero  and  updated  at  the  end  of  each  pass  through  the  routine  as 
well.  These  two  variables  are  packed  as  shown  in  Table  IV-8  in  integer  vari¬ 
ables  IFLOD^.  and  IPAS  ,  where  1=1  ,  KFT  ,  the  number  of  flood/drv  cells. 

Table  IV-8 

Flood/Dry  Logic  Control  Variables 


I 


IFLOD,  =  IFL  +  (109)  IQU  +  ( 106)  IQV 

cell  face”^ 
counter 


v  cell  face 
counter 


IFL  =  1000  *  N  +  M 

Crid  cell  indices 


v  face  status 
indicator 


ipast  =  f.o1) 


IPS  +  JP£ 
tu  face  status 


Indicator 


213.  Three  distinct  cell  face  states  are  distinguished: 

IPS  =  0  Cell  face  fully  open  or  initially  closed 

IPS  -  1  Cell  face  partially  open  one  side 


IPS 


Cell  face  partially  open  both  sides 
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ISP  =  3  Cell  face  closed  after  being  fully  open 
The  u-face  opening  logic ' proceeds  as  follows: 

1.  Do  not  attempt  to  open  a  cell-face  previously  closed,  unless 

k+1 

it  has  been  closed  for  KS1  time  steps.  Define  T1  =  - 
,  T9  k+1  n,m 

n  ,m+l 

2.  For  T1  and  T2  >  DFLUE  ,  IPS  =  2  ,  perform  the  procedures 
shown  in  Table  IV-9. 

3.  For  either  T1  or  T2  >  DFLUE  ,  IPS  =  1  ,  perform  the  fol¬ 
lowing  procedure;  illustrated  for  the  case  TL  >  DFLUE  . 


£ 


1/2  9 

Q  =  minjCd(Tl  -  DFLUE)  ,  (uj )  (Ac^ ) “ (i^ )  L | T1 


T2|/^Ac((u1)m  +  Cu^J] 


An,  =  -Qt/ (u, )  Act, 
i  i  m  l 


An2  =  QT/Cu^^jAaj 


k+1  k+1 

=  n  +  An 
n,m  n , m  L 


k+1  k+1 

,  ,  =  n  +  An9 

n , m+ 1  n , m+ 1  2 


Analogous  procedures  hold  for  the  case  T2  >  DFLUE  . 

4.  If  T1  and  T2  <  DFLUE  ,  consider  the  next  flood  cell  face. 

215.  Cell  face  conditions  are  then  examined  for  closing.  The  u-face 
:Iosing  logic  proceeds  as  follows: 


1 


With  TI  and  T2  defined  as  previously,  if  T1  and 
T2  >  DFLU-EPSD/2  ,  consider  the  next  flood  cell  face. 

If  Tl  or  T2  or  both  are  less  than  DFLU  -  EPSD/2  , 
consider  the  cell  face. 


ICU 


II  = 50 

t 

UPn  "  0 
*. 

IQU  =  1 

IPS  =  3 

n,  m 

n  ,m 

1  58 
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T 1  -  OFLuE  :•  0 
T2  -  OFluE  • Q 


sp  -  oeta  -  n 
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REFERENCE 
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SO  --  OETA  -  T2  <0 


•  Jf-x 

T2 f  TDA  -  12  IT!  *  T2I  -  DFLUE  I AVERAGE  FLOOD  DEPTH) 
\OETAI^  Jj)FLUE 


DETA  -  Tl  =  SP  <  0  I 

r 

i 

■  "JfOETA'^l  SO  =  DETA  -  T2  >  0 

T’  J 

TDA 

n 

T2 

f  OF L UE  MODEL  DATUM 

-L  REFERENCE 
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DETA  =  d 


dlTI  +  d2T2 
d,  +  d. 


S?  =  d  -  Tl  <  0  ,  SO  =  d  -  T2  >  0 

n 


(Tl  +  T2) 

IDA  =  —  ,  —  -  DFLUE 


j- 


0  =  nin  ,  (TDA) 
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3/ 


’  ■yl)m<-‘i0‘l)  U 1  ^m+1 1 
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If  T12  =  0  consider  Che  next  flood  cell  face 


1  m+1 


(Continued) 
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Table  IV-9  (Continued) 


'1. 


Til  •  0 


d..  =  -Q t / d <  0 

If  d_  <  SQ  ,  then  use  SP  and  SQ 


d  =  -o-r/d,  <  0 
'l  11 


as  computed  above  for  corrections. 
If  d_  >  SQ  ,  compute 


If  d^  SP  ,  then  use  SP  and  SQ 
"  I 

as  computed  above  for  corrections. 

If  d^  >  SP  ,  compute 
"l 


d-  =  -d„  *  d,/d.  . 

1  ■“*  W  1. 


d  =  -d  *  d, /d^  . 

n„  n .  II 


1 


?et  SP  ^  d 


SO  -  d 


IF  (SP  +  Tl)  consider  next  flood  cell  face;  i.e.,  do  not  allow  negative 

water  depth  for  cell  II. 


IF  (SQ  +  T2)  <_  consider  next  flood  cell  face;  i.e.,  do  not  allow  negative 


water  depth  for  cell  MM1, 


Else,  SEM 

=  SP 

+  semtt 

II 

n 

>  Note 

corrections 

QFM 

5  “MM1 

=  S0 

+  SEV 

IPS 

=  I 

(Cell  starting 

to  open) 

ext  compute  cell  face  velocitv: 


,y  =  _ ±2i 

•  (u  ) 


_k+l  +  „k 
i  i,  \  n  , m+ 1  n  , m+ 1 


1 i .  \  n , m+ 1 

m+I/2  1 


_  k+I  __k 

n  ,  m  n  ,  m 


DPTH  =  max 


0. 5 ("k  .  -  h  +  -  h  | ,  EPSD2 

\  n , m+ l  n,m+I  n,m  n ,  m 
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Table  IV-9  (Concluded) 


Q'  =  DPTH  *  (U  +  Q') 

If  i Q ' i  >  Q  »  then  Q  =  sign  (a,  Q') 
t'Pn  =  Q/DPTH 


If 


(SP  +  Tl)  <  H  +  EPSD  \ 

l  Consider  next  flood  cell  face 

(SQ  +  T2)  <  Hmm1  +  EPSD 


Else,  set  cell  to  fully  open  (normal)  status 


ICU  =  60 


IQU  =  1  ,  IPS  =  0 


216.  As  in  the  revised  technique  of  Butler  and  Prater  [37],  all  water 
surface  elevation  adjustments  are  accumulated.  Therefore,  the  cell  face  com¬ 
putations  are  order  dependent.  Thus  in  alternate  time  steps  u  faces  are 
considered  before  v  faces  and  vice  versa.  In  addition,  the  order  of  the 
potential  flood/drv  cell  list  is  reversed  each  time  step.  The  drainage  con¬ 
cept  has  been  eliminated  in  the  interest  of  computational  economv.  In  addi¬ 
tion,  only  cell  face  closing  are  considered  for  both  u  and  v  faces  each 
time  step  as  indicated  as  follows: 

Odd  Time  Step  Number: 

u  face  opening 
u  face  closing 
v  face  closing 
Even  Time  Step  Number: 

v  face  opening 
v  face  closing 
u  face  closing 


217.  If  during  the  computations  in  Subroutine  MOTN  a  negative  water 
depth  is  encountered  in  a  flood/dry  cell,  this  cell  is  immediately  removed 
from  the  computation  by  setting  its  Chezv  C  and  all  three  time  levels  of  its 
water  cell  face  velocities  to  zero.  In  addition,  all  three  time  levels  of 
water  surface  elevation  are  set  equal  to  the  cell  bottom  elevation.  The 
amount  of  water  at  time  level  k-1  is  accumulated  and  distributed  equally 
over  all  wet  cells  not  on  the  boundary  as  an  adjustment  on  most  recentlv  com¬ 
puted  water  surface  elevation.  A  critical  Courant  number  Is  computed  and  out 
put  to  indicate  potential  difficulty.  Within  Subroutine  UPDATE,  it  is  possi¬ 
ble  to  add  EPS;)!  feet  of  water  to  non-flood  cells  which  develop  negative 
water  depth.  If  this  occurs  for  several  cells,  then  DSCOl'R  needs  to  be 
further  decreased  in  order  to  include  these  cells  as  flood  cells. 


development  of  the  numerical  imple- 
-er.tatlort  of  the  short  wave  equation 

218.  The  short  wave  equations  (2.  1.6-2.  1.7  of  Part  III  and  2.2.1  and 
2.2.2  of  Part  III)  require  the  specification  of  a  fetch  length  and  average 
water  depth.  In  hurricane  events  which  pass  over  Lake  Okeechobee,  the  lake 
windfields  are  highlv  curved.  The  following  approach  was  developed  to  develo 
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an  effective  fetch  in  such  cases.  In  addition,  the  dynamic  Land/water  bound¬ 
ary  must  be  taken  into  account  in  order  to  develop  fetch  lengths.  It  is  as¬ 
sumed  that  the  tree  islands  or  other  subgrid  barriers  do  not  influence  the 
windfield  sufficiently  to  be  taken  into  account  in  the  fetch  determination 
procedure  as  outlined  below. 

a.  Compute  ZTRY(J.I)  array  for  J  =  L  ,  NLAX  ,  and  1=1, 
MLAX  over  the  uniform  hurricane  auxiliary  grid.  If 
ZTRY(J,I)  is  greater  than  or  equal  to  0.5  cell  (J,I)  is 
considered  dry. 

b.  Define  the  following  eight  ordered  direction  pairs: 

Octant  No: _ 1 _ 2 _ 3 _ 4 _ 5 _ 6 _ 7 _ 8 

IXKS  1  1  0  -L  -1  -1  0  1 

1YKS  0  1  1  1  0  -1  -1  -I 

c.  ■  Wave  calculations  are  performed  along  all  hydrodynamic  grid 
faces  along  the  levee.  The  effective  depth  used  in  the  wave 
equations  is  taken  as  the  water  depth  at  that  time  in  the  cell 
of  the  WIFM  grid.  Since  grid  cells  are  on  the  order  of 
0.35  square  miles  in  size,  and  the  wave  height  may  be  depth 
limited,  this  procedure  is  considered  acceptable.  It  may  be 
desirable  to  consider  the  effect  of  alternative  depth 
considerations . 

d.  In  order  to  integrate  the  fetch  calculation  procedures  with 
the  ZTRY  array  logic,  it  is  necessary  to  transform  the  loca¬ 
tions  of  levee  subgrid  barriers  on  the  WIFM  hydrodynamic  grid 
to  the  uniform  hurricane  auxiliary  grid.  Consider  a  levee 
section  transformed  to  cell  (J*,  I*)  on  this  uniform  grid. 

The  fetch,  F  ,  is  computed  as  follows: 

(1)  Set  F  =  0  . 

(2)  Determine  the  direction  octant  number  from  which  the  wind 
is  blowing. 

(3)  Determine  the  appropriate  direction  order  pair  IXKS*  , 
IYKS*  . 

(4)  Compute  a  fetch  length  scale,  ASp  ,  where 
AS  =•  yj( IXKS*)2  +  (IYKS*)2  AS 

and  AS  is  the  uniform  grid  spacing. 

(5)  F  =  IAS_  ,  where  I  is  such  that 
F 

ZTRY ( J*  +  ilYKS*  ,  I*  +  ilXKS*)  <  0.5 
for  1  <  i  <  I  and 


ZTRY  ( J*  +  (I  +  1 )  IYKS*  ,!*+(!+  DIXKS*)  >  0.5 


219.  The  fetch  is  considered  zero  if  the  wind  magnitude  is  zero.  The 
implementation  of  the  runup  and  overtopping  computational  procedures  were  es¬ 
sential  v  performed  by  Mr.  David  L.  Leenknecht,  while  he  was  at  the  New 
Orleans  District.  Several  curves  were  computerized  as  shown  in  the 
Table  1V-10.  The  wave  computations  are  performed  every  KS1  time  steps.  In 
present  applications,  KS1  =  L2  corresponding  to  a  1  hr  interval  for  a  hydro- 
dynamic  time  step  of  300  seconds.  Wave  overtopping  volumes  are  not  written  on 
the  disc  transfer  file  for  input  to  the  Link-Node  drainage  model  in  the  pres¬ 
ent  model  version.  However,  these  volumes  are  computed  and  are  written  on  the 
output  file. 

Table  IV-10 

Wave  Computation  Routines 


Subroutine  or  Function  Name 


Function  Eq.  8 
Function  Fig.  75 
Function  Eq.  7 
Function  Fig.  50 
Function  Eq .  79 
Subroutine  RL'NTP 
Function  Fig.  74 
Function  Fig.  75 


Purpose 


Eq.  8  Fig.  20 
Figure  7-5 

Overtopping  equation 
Figure  50 

Overtopping  equation 
Figure  6,  Pg.  25 
Figure  7-4 
Figure  5 


Technical  Source* 


TP-78-2  Stoa  [32] 
SPM  [27] 

J.  R.  Weggel  [33] 
TP-78-2  Stoa  [32] 
J .  R.  Weggel  [  33 ] 
TP-78-2  Stoa  [32] 
SPM  [27] 

TP-78-2  Stoa  [32] 


220.  The  major  computational  routines  comprising  the  model  and  their 
purpose  are  as  follows: 

a.  Subroutine  WAVES - Shallow  water  model  driving  routine 

b.  Subroutine  SWAVEL - Shallow  water  wave  height  and  period 

(i)  CETN  -  I  -  6  [28] 

(ii)  CW-167  Curves  [31] 

c.  Subroutine  WAVED - Newton-Raphson  determination  of  the  shallow 

water  wave  length 


d.  bubrouting  DEEP - Determination  of  the  shoaling  factor  and 

deepwater  wave  height  and  length 

dll.  The  incorporat ion  of  short  wave  computations  directly  within  the 
long-wave  model  represents  an  extension  over  previously  employed  methods. 
Within  the  federally  supported  flood  insurance  studies,  short  wave  maximum 
heights  were  computed  ipso  facto  after  the  long  wave  peak  surge  values  were 
known.  The  fetch  determination  procedures  as  well  as  the  wind  magnitude  were 
not  as  accurate  as  the  methods  developed  here  for  Lake  Okeechobee. 
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PART  V:  SEICHE  CALIBRATION  AND  VERIFICATION  OF  THE  LONG  WAVE  MODEL 


Global  grid  development 

222.  In  order  to  study  long-wave  phenomena  on  Lake  Okeechobee,  it  is 
necessary  to  construct  a  grid  on  which  the  hydrodynamic  computations  are  car¬ 
ried  out.  The  resolution  of  the  grid  determines  the  time  step  and  the  number 
of  cells  over  which  computations  are  performed.  In  order  to  forecast  hurri¬ 
cane  phenomena  for  various  possible  storm  intensities  and  tracks,  it  is  neces¬ 
sary  to  construct  a  grid  such  that  simulations  on  the  order  of  24-36  hours  not 
be  cost  prohibitive.  For  this  purpose,  a  64  *  48  (3072  cell)  grid  was  con¬ 
structed  as  shown  in  Figure  V-l.  The  grid  spacing  is  uniform  in  the  horizon¬ 
tal  and  variable  in  the  vertical  with  larger  spacings  in  the  center  of  the 
lake,  where  less  resolution  is  required.  Each  direction  is  mapped  indepen¬ 
dently  using  a  time  sharing  code  (Program  MAPIT) ;  the  results  are  shown  in 
Table  A-l  (Appendix  A). 

223.  The  maximum  spatial  resolution  is  3333.3  ft.  For  a  gravity  wave 
speed  of  approximately  20  fps ,  an  explicit  time  step  of  180  seconds  is  more 
than  adequate.  In  the  alternating  direction  implicit  scheme  employed,  a  time 
step  of  360  seconds  is  used  corresponding  to  a  Courant  number  of  two. 

224.  The  geometrical  and  topographic  characteristics  of  the  Lake 
Okeechobee  system  are  resolved  on  the  global  grid  developed  above. 

Jacksonville  District  (SAJ)  personnel  digitized  the  topography  as  shown  in 
Table  A- 2  employing  nautical  chart  L1428,  18th  Edition,  Aug.  23/80.  All  land 
elevations  are  with  respect  to  mean  sea  level  and  are  shown  as  positive  num¬ 
bers.  At  the  outermost  boundary  of  the  grid  an  arbitrary  elevation  of  59  feet 
was  specified.  This  elevation  is  such  that  no  water  elevation  will  exceed  it 
even  in  worst  case  hurricane  simulations;  thus,  all  water  will  be  contained  on 
the  computational  grid.  Water  depths  are  specified  with  respect  to  a  datum 
plane  located  12.5  ft  above  mean  sea  level  and  are  preceded  by  a  minus  sign. 
The  model  datum  plane  is  taken  as  the  depth  datum  plane.  The  levee  system  and 
tree  islands  are  specified  as  barriers  occupying  a  complete  side  of  a  grid 
cell.  Their  geometrical  and  physical  characteristics  were  provided  by  SAJ  and 
are  presented  in  Tables  A-3  and  A-4 ,  respectively.  The  general  configuration 
of  the  levees  and  Tree  Islands  are  shown  in  Figure  V-2.  Gage  locations  where 
water  surface  elevation  data  are  available  are  also  as  indicated  in  Figure  V-2 


Seiche  Phenomena 

22 3.  Consider  a  spatially  uniform  wind  acting  over  the  lake  for  a  given 
period  of  time.  The  wind  exerts  a  surface  stress  on  the  water  in  the  direc¬ 
tion  of  the  wind.  During  the  seiche  generation  phase,  the  water  sets  up  on 
the  downwind  side  of  the  lake  and  sets  down  in  the  upwind  side  of  the  lake. 

If  the  wind  remains  constant  for  a  sufficiently  long  period  of  time,  a  so- 
called  static  wind  tide  is  established.  The  term  "static"  is  used  to  denote  a 
condition  in  which  all  transient  water  level  fluctuations  due  to  the  time 
variation  and  start-up  of  the  wind  have  been  damped  out  by  the  lake  frictional 
mechanisms . 

226.  During  the  seiche  decay  phase,  the  wind  approaches  zero  and  the 
static  water  level  surface,  oscillates  until  the  lake  frictional  mechanism 
removesall  the  potential  energy  associated  with  the  initial  water  surface 
displacement . 

227.  This  hypothetical  seiche  phenomena  is  sketched  in  Figure  V-3.  In 
actuality,  the  wind  in  the  seiche  generation  phase  may  not  be  spatially  uni¬ 
form  and  the  wind  magnitude  may  not  remain  constant  for  a  sufficiently  long 
period  of  time  to  develop  a  static  wind  tide.  During  the  seiche  decay  phase, 
the  wind  magnitude  will  usually  be  low  but  not  completely  zero. 

228.  Consider  the  physical  mechanisms  occurring  in  each  phase  of  the 
seiche.  In  the  generation  phase,  the  wind  exerts  a  surface  stress  on  the 
water.  The  bottom  stress  is  also  in  the  wind  direction;  i.e.,  in  the  same 
direction  as  the  surface  stress  due  to  the  return  flow  (flow  reversal)  at  the 
bottom.  In  a  two-dimensional  vertically  integrated  model,  the  return  flow 
phenomena  cannot  be  simulated.  The  velocity  is  integrated  over  the  vertical 
to  produce  an  average  velocity.  The  bottom  stress  occurs  in  the  direction 
opposite  to  the  averaged  velocity  vector.  During  the  static  tide  phase,  the 
surface  wind  force  is  balanced  by  the  pressure  gradient  forces  and  the  veloc¬ 
ity  in  the  water  bodv  approaches  zero;  therefore,  there  are  no  bottom 
stresses.  During  the  decay  phase,  the  wind  surface  stress  approaches  zero; 
the  pressure  gradient  force  supplies  the  energy  to  move  the  water;  and,  this 
energy  is  dissipated  by  the  bottom  stress. 

229.  In  order  to  calibrate  the  bottom  stress,  it  is  necessarv  to  simu¬ 
late  the  seiche  decay  phase,  since  it  represents  the  only  surface  stress 
acting  on  the  water  body.  In  order  to  calibrate  the  wind  surface  stress,  it 
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is  necessarv  to  simulate  the  static  tide  case,  since  then  the  only  surface 
stress  applied  to  the  water  body  is  due  to  the  wind.  Once  both  the  bottom 
friction  and  surface  wind  stresses  have  been  calibrated,  it  should  be  possible 
theoret ical ly  to  simulate  the  complete  seiche  (generation  static  tide,  and 
decay  phase).  However,  in  practice,  the  difference  between  the  static  tide 
level  and  the  initial  water  surface  at  a  given  gaging  station  is  approximately 
.5  feet.  In  addition,  the  initial  water  surface  elevations  at  any  given  time 
are  not  equal  at  all  gaging  stations.  The  vertical  control  was  accomplished 
through  first  order  leveling  and  is  probably  accurate  to  .1  feet.  The 
spec  if icat ion  of  the  initial  water  surface  is  further  complicated  by  the 
availability  of  data  at  only  6-7  stations.  The  above  constraints  make  it 
extremely  difficult  to  simulate  the  complete  seiche.  As  a  result,  the  work 
presented  here  focuses  on  the  seiche  decay  phase  and  the  calibration  and  veri¬ 
fication  of  the  bottom  friction  mechanics.  The  August  12-13,  1968  seiche  was 
used  to  calibrate  the  bottom  friction  mechanics,  while  the  August  20-21,  1964 
seiche  was  used  to  verify  the  calibrated  friction.  These  two  periods  were 
considered  by  Whitaker  and  Reid  [40]  in  their  numerical  study  of  Lake 
Okeechobee  and  were  selected  so  that  the  results  of  the  present  study  could  he 
compared  with  this  previous  work. 

August  12-13,  1968  seiche  calibration 

230.  Wind  and  stage  data  were  obtained  from  the  Jacksonville  District 
it  the  locations  shown  in  Table  A-6 .  Measured  water  surface  elevations  at 
MGS-2  and  HGS-6  are  shown  in  Figure  V-4 .  The  period  from  2000  HRS  on  August 
August  12  through  1200  HRS  on  August  13  was  simulated.  This  period  corre¬ 
sponded  to  the  seiche  decay  phase.  The  wind  data  are  presented  in  Table  A-7. 
As  noted  from  Table  A-7,  all  wind  magnitudes  are  less  than  3  mph  except  for  a 
tew  readings  at  Port  Mayaca.  Water  surface  elevation  characteristics  are 
given  in  Table  A-8 .  Average  station  water  surface  elevations  are  within  a 

.1  feet  except  at  LS-12,  which  appears  to  be  .1-.2  feet  high.  Recorded  lake 
station  elevations  correspond  extremely  well  with  the  outside  gage  readings  at 
til  stations  over  the  measured  period.  Gage  timing  information  is  shown  in 
Table  A-9.  All  recorded  timings  are  within  1-2  hours  over  a  period  of 
ipproximate ly  one  month. 

231.  In  order  to  account  for  vegetation  effects,  the  following  effective 
bott  5m  I r it  ion  mechanism  is  calibrated: 


/  n. 
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where 


n*  =  effective  Manning's  n  roughness 
nd  =  Manning’s  n  dependent  upon  water  depth  with  respect  to  model 
datum,  d 


cd 


I  =  Canopy  coefficient  1  with  respect  to  water  depth  d 


cd9  =  Canopy  coefficient  2  with  respect  to  water  depth  d 
dean  =  Canopy  friction  depth 

The  relations  given  in  Table  A-10  were  calibrated.  The  factor 

(1.  +  cd^e  d/cd2  j  is  evaquated  in  Table  A-ll  for  various  values  of  depth  less 

than  or  equal  to  the  canopy  friction  depth. 

232.  In  order  to  illustrate  the  friction  specification,  consider  water 
cell  (33,  33).  From  Table  A-2,  the  water  depth  with  respect  to  model  datum  is 
eight  feet.  Consulting  the  friction  versus  depth  relationship  table,  we  note 


eight  feet  is  in  depth  range  4,  and  n.  =  .018.  From  the  next  table,  the 


canopy  multiplication  factor  for  8  feet  equals  1.5413.  Therefore,  the 
effective  roughness  n*  for  cell  (33,  33),  equals  (1.5413)  x  ,018  =  .0277. 
This  value  is  used  throughout  the  course  of  the  simulation,  due  to  the  fact 
the  water  depth  remains  essentially  constant  in  seiche  computations. 

233.  In  order  to  specify  initial  conditions,  the  following  approach  was 
used.  All  velocities  were  assumed  equal  to  zero.  The  initial  water  surface 
as  shown  in  Table  A-12  was  developed  by  employing  Subrouting  TILT  as 
documented  in  Appendix  D.  For  the  August  12-13,  1968  Seiche,  a  plane  was  used 
to  specify  the  initial  water  surface.  The  plane  passed  through  the  Clewiston 
HGS-2  recorded  level  (15.20  ft)  and  through  the  Okeechobee  HGS-6  recorded 
level  (l-*.52  ft)  at  2000  HRS  on  12  August.  No  lateral  tilt  was  used.  The 
measured  data  at  half  hour  intervals  are  shown  in  Table  A-13  from  2000  HRS  12 


measured  levels  in  Plates  B-l  through  B-7  (Appendix  B). 

In  interpreting,  the  ordinate  0.00  NGVD  corresponds  to  15.0  ft  MSL.  In 
the  top  half  of  the  Plates,  the  simulated  levels  not  including  the  tree 
islands  are  shown.  Simulation  results  including  the  tree  islands  are 


snovn  in  the  lower  half  ot  the  plates  for  the  sake  t  comour ison .  At  most 
stations,  the  tree  islands  appear  to  exert  a  vt-rv  '.mall  int  lucnce  on  the  water 
levels  during  the  seiche.  Simulated  water  level-,  correspond  well  with  mea¬ 
sured  levels  at  all  stations  except  LS-L2.  As  noted  from  Table  A-S  ,  the  mean 
at  LS-12  is  . 1 - . 2  ft  higher  than  at  all  other  stations.  This  suggests  that 
there  may  be  a  geodetic  leveling  problem  at  I.S-12  of  . 1  - . 2  feet;  i.e.,  it  mav 
be  appropriate  to  subtract  .1-.2  feet  from  the  recorded  level  at  iS-12.  If 
this  is  done,  the  simulated  levels  would  match  the  adjusted  levels  verv 
closely  at  LS-12. 

August  20-21,  1964  seiche  verificat ion 

235.  Wind  and  stage  data  were  obtained  from  the  Jacksonville  District 
at  the  locations  shown  in  Table  A-13.  Measured  water  surface  levels  at  HGS-2 
ind  HGS-6  are  shown  in  Figure  V-5  in  order  to  illustrate  the  seiche  phenomena. 
The  period  from  1730  HRS  on  August  20  through  0900  HRS  or.  August  21  was  simu¬ 
lated.  As  may  be  determined  from  the  wind  data  presented  in  Table  A-12,  this 
corresponded  to  the  seiche  decay  phase.  All  wind  magnitudes  are  less  than 
o  mph  throughout  most  of  the  period  at  all  gages  extent  for  a  few  hours  at 
Port  Mayaca.  Water  surface  elevation  character ist ics  are  given  in  Table  A-13. 
Average  station  water  surface  elevations  are  within  .1  foot  except  at  HGS-2 
and  LS-ln.  LS-16  appears  to  be  .03-. I  feet  high.  Recorded  lake  station 
elevations  correspond  extremely  well  with  the  outside  gage  readings  at  all 
stations  over  the  measured  period.  Gage  timing  information  is  shown  in 
Table  A-lo.  HGS-2  and  HGS-n  were  reset  prior  to  the  simulated  seiche  and 
should  be  extremely  accurate  timewise.  LS-1-*  and  In  exhibited  fairly  sub¬ 
stantial  '  4 . 3-3 . 3  hrs)  timing  errors  in  contrast  to  the  Part  Mayaca  .3  hr  error 
veer  i  vie  month  period. 

2  3n .  Th<-  et  feet  ive  bottom  friction  mechanism  in  liquation  (2.1.1)  and 
the  depth  versus  friiti.ci  range  re  lat  ion  in  Table  A-10,  and  the  canopy  coef- 
; i .  ient  s  in  Tibie  A- 1 1 ,  were  use a  in  the  simulation.  The  initial  water  surface 
was  is  i  W.  eve  J  pari:  1  i  is  ut  lined  in  Appendix  B  with  no  lateral 

til:  p  i  g  :  ,r  .gn  •  nt  s  »h  v.  :::  i  r  i  e  A- 1  7  .  The  measured  data  at  half 

.  r  .  e  :  o  vu  .  1  r  .•  ;r  m  1  '  >  :  HRS  20  August  to  u900  HRS 

\  .  .  •  .  ;t...  i  •  •  :  -p  ir,  :  versus  measured  levels  in 


1.0  HGVD  corresponds 
1  ■.  •  e  i  levels  not 


including  the  tree  islands  are  shown.  Simulation  results  including  the  tree 
islands  are  shown  in  the  lower  half  of  the  plates  for  the  sake  of  comparison. 

At  the  majority  of  the  stations,  the  tree  island  appear  to  exert  a  very  small 
influence  on  the  shape  of  hydrographs.  Simulated  water  levels  correspond  well 
with  measured  levels  at  all  stations  except  at  LS-lb.  As  noted  from 
Table  A-15,  a  geodetic  leveling  problem  at  LS-16  on  the  order  of  .1  feet  may  be 
indicated.  If  .1  feet  is  subtracted  from  the  recorded  levels  at  LS-lb,  the 
simulated  levels  would  more  nearlv  correspond  to  these  adjusted  levels. 


PART  VI;  AUGUST  1949  HURRICANE 


General  Ftorm  Characteristics 

237.  The  August  26-27,  1949,  hurricane  was  used  to  verify  the  previ¬ 
ously  calibrated  bottom  friction  relationships  developed  in  the  seiche  work 
presented  in  Part  V.  Two  separate  hurricane  models  were  used  to  develop 

I 0-meter  10-minute  average  sustained  overwater  windspeeds  for  use  in  the 
hydrodvnamic  sub-model.  The  Powell  wind  stress  formula  [23]  was  used  to  com¬ 
pute  surface  wind  stress  in  the  hydrodynamic  sub-model  from  the  10-meter 
ln-minute  average  windfields.  The  overwater  windfields  were  reduced  over  Lake 
Okeechobee  dynamically  based  upon  land/water  boundaries  developed  in  the 
hydrodynamic  simulation. 

238.  Initially  results  are  presented  for  the  two  hurricane  sub-models 
in  terms  or  10-meter  10-minute  average  overwater  winds.  The  August  26-27, 

'.Q4Q  hurricane  with  a  maximum  central  pressure  intensity  of  1.8  in  Hg  repre¬ 
sents  the  most  severe  storm  to  have  influenced  Lake  Okeechobee  for  which  mea¬ 
sured  «ater  level,  wind  and  pressure  data  are  available.  Therefore,  in  a 
studv  of  wind  generated  water  levels  on  the  lake,  this  storm  must  be  con¬ 
s’.  iered.  The  track  of  the  storm  is  shown  in  Figure  VI-1.  As  may  be  observed 
r'-e  storm  passed  directly  over  the  northern  section  of  the  lake.  Due  to  the 

■.r.ter  i'ckwise  circulation  around  the  low  pressure  center  and  the  track, 

.•  i  n  i  direct  i  'ns  shift  draraat  leal  ly  f  ram  out  of  the  North  to  out  of  the  South 
in*  the  :  curse  of  the  storm. 

:  me  sub — ode  1 

-  it  r.  results  1  CPH  and  PH!.' 

* .  In  both  hurricane  simulations,  the  position  of  the  storm  is  as 
'  e  -  :n  "  »h  *.  e  » - .  '  in  hours  alter  the  initial  position  corresponding  to 

-eenwich  mean  time  on  August  26 ,  1R49.  Thirtv-one  hours  oi  wind- 

•  i  e .  n:  rmatim  is  '  vrputed  emp  loving  five  snapshots.  Snapshot  data 

e~r '  e  :  in  the  !”:!  h.urricane  simulation  are  presented  In  Table  —  2 1  .  cnap- 
-■ .  t  v  1  r  *  f  1  e  1  i  results  are  presented  In  i’latvs  3-13  through  l.  -  1  7  in  graphical 

•  ’T — .  ‘he  entire  spatial  windfleld  Is  presented  and  centered  with  respect  to 
:ue  eve  in  these  figures.  'hserve  the  inflow  angle  Is  constant  at  each  radial 
‘.'.stance. 


Figure  VI-1.  August  1949  Hurricane  Track 


240.  The  isovel  pacterns  are  shown  in  Places  B  —  1 8  through  B-22  for 
each  corresponding  snapshot.  The  wind  structure  only  on  the  first  two  inner 
nose  grids  of  the  five  level  nest  is  shown  in  these  figures.  Single  digit 
numbers  represent  isovel  contours  at  a  factor  of  ten  in  knots.  Snapshot  data 
employed  in  the  PBL  hurricane  simulation  are  shown  in  Table  A-21.  Note  the 
data  correspond  to  those  employed  in  the  SPH  simulation  with  the  exception  of 
the  steering  flow.  Steering  flow  is  not  used  directly  in  the  SPH  simulation. 
In  the  PBL  simulation,  a  shearing  flow  of  9  meters/sec  is  specified  in  the 
direction  of  storm  motion.  This  represents  a  typical  value  for  storms  in  this 
region.  Snapshot  windfield  results  are  presented  in  Plates  B-23  through  B — 2 7 
in  graphical  form  in  the  same  format  as  for  the  corresponding  SPH  results. 
Observe  the  inflow  angle  is  computed  dynamically  in  the  PBL  simulation  and  is 
no  longer  constant  at  each  radial  distance.  Isovel  patterns  are  given  in 

SPH  corresponding  format  in  Plates  B-28  through  B-32  for  the  PBL  simulation. 
The  resulting  structure  is  quite  different  than  the  SPH  result.  Maximum 
gradient  winds  (kts)  and  maximum  and  PBL  simulations  in  Table  A-22  for  each 
snapshot.  PBL  simulations  results  are  generally  3-8  knots  less  than  SPH 
maximum  wind  results. 

241.  Since  the  land/water  boundary  influences  the  reduction  factor 
array  to  adjust  the  10-m  10-minute  average  windfield,  which  in  turn  drives  the 
hydrodynamics  determining  the  land/water  boundary,  in  effect  a  dynamic  feed¬ 
back  loop  is  being  simulated.  It  is  necessary  then  to  consider  the  hydro- 
dynamic/windf ield  computations  as  a  coupled  system. 

2^2.  Results  for  the  windfield  solution  and  pressure  field  at  gige 
locations  shown  in  Figure  VI-2  are  presented  in  Plates  B  —  3 3  through  3--n  for 
the  PBL  simulation  and  in  Plates  B-47  through  B-60  for  the  SPH  simulation. 

The  double  hat  wlndspeed  distributions  computed  at  HGS-b,  Port  Mayaca,  I.S-1J, 
ind  LS-19  ire  in  agreement  with  observed  profiles  shown  in  CW-167  [31].  The 
central  pressure  distributions  are  the  same  in  both  SPH  and  PBL  simulations 
and  correspond  to  the  Schloemer  profile.  The  spatial  structure  of  the  wind- 
field  in  the  vicinity  of  Lake  Okeechobee  is  shown  in  Plates  B-61  through  B-bH 
for  the  PBL  and  in  Plates  B-69  through  B-76  for  the  SPH  simulation.  Observe 
how  the  windfield  increases  overwater.  At  simulation  time  17  hours 
corresponding  to  2300  HRS  EST  on  26  August,  the  placement  of  the  eve  in  the 
nor theastermost  corner  of  the  lake  is  excellent  agreement  with  placement 
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data  shown  in  C*-lb7  [31].  Note  the  general  wind  direction  and  structure  ire 

similar  except  near  the  eve  of  the  storm.  HourLy  predicted  windspeed  and 

direction  data  for  both  PBL  and  SPH  simulations  are  compared  versus  measured 

data  in  Table  A-23.  For  the  first  3-9  hours  PBL  and  SPH  winds  are  in  close 

lgreement.  However  from  hour  Id  through  22  SPH  winds  are  as  much  as  JO-13  kt.-. 

stronger  at  some  gage  locations.  Near  the  end  of  the  storm,  results  are  in 

closer  agreement.  This  suggests  that  in  the  near  eye  region  the  general 

structure  of  the  windtields  are  considerably  different  but  .ire  essential Lv  the 

same  in  the  storm  regions  greater  than  two  to  three  radii  to  maximum  winds. 

A-tual  hydrodynamic  results  are  presented  in  turn  for  each  hurricane 

s imu 1  at  ion . 

PBL  hvdrodynamic 

sun-model  simulat ion  resu  Its 

J 1 3 .  A  period  of  J8  hours  starting  at  ObOO  HRS  EST  on  2b  August  19.9 
was  simulated  using  a  time  step  of  130  seconds.  Convective  acceleration  terms 
were  onsidered  in  the  middle  of  the  lake.  The  previously  presented  PBL  wind- 
:  ield  md  Schloemer  pressure  t'ieLd  were  used  to  derive  the  hydrodynamics.  No 
i:u  lows  ’  r  rainfall  were  considered.  The  levee  coni  igurat  ion  for  the  August 
burr:  ar.e  is  presented  in  Figure  VI-3.  Simulated  water  level  histories 

■  re  it"  :  with  hi  served  water  level  histories  in  Plates  B  —  7 7  through  B-8  3. 

1  tt  ions  m  Plates  3-77  through  B-3  5  are  with  respect  to  an  initial  like 

...  1  !  1  . .  ;  :  t  ms  1  (plot  zero)  .  Simulated  levels  are  in  reasonable  good 

tgr-  e-.e-it  with  observed  Levels  except  at  HGS-3,  HCS-b,  and  LS-lb.  The  simu- 
!  .  t  HGS-*>  i  Okeechobee)  is  approximately  4  feet  below  the  observed 

•  o  .  ••  t iming  t  the  simulated  peaks  and  troughs  in  water  level  records  is 

■:-  r  i.  .  ;u  1  «.-»«•  agreement  with  the  observed  data.  This  fact  substantiates 
•  ,-rii  -n  track  and  windtield  structure  produced  bv  the  PBL  simulation, 

it.-.  1  .  •  :  .r  rent  histories  (vertical  lv  integrated)  are  shown  in  Plates  B-Sb 

■  or  .ci  j-->2 .  ib  acral  current  structure  o.er  the  entire  lake  is  represented 

.  ;  ■  ’ !  i ■  •  ->  i  through  3- 100.  Maximum  current  magnitude  is  on  the  order  of 
j.  .  iim.il  i  ted  significant  wave  heights  and  periods  are  compared  with 

lit  i  i:i  1  is!--  A-J  .  .  Maximum  simulated  significant  wave  heights  are  on  the 

rier  it  1  teet  with  corresponding  periods  of  approximate  lv  ',.3  econds. 
Maximum  -.till  water  levels  are  presented  in  Table  A-25.  Boxed  cells 
represent  t 1. coded  areas.  Areal  extent  of  simulated  flooded  areas  is  in 
lgreement  with  d it  a  presented  in  CW-lb7  [31].  Minimum  stilled  water  depths 


are  shown  in  TabLe  A-’b.  Boxed  cells  represent  areas  initially  wet, 

which  dried  sometime  luring  the  simulation.  A  minimum  depth  of  1  foot  or 

Less  was  used  as  the  demarcation.  The  areal  extent  of  exposed  areas  shown  is 

in  general  agreement  with  data  presented  in  CW-167  [31]. 

SPH  hydrodynamic 

sub-model  simulation  results 

The  same  28  hour  period  was  simulated  employing  a  time  step  of 
130  seconds.  The  previously  presented  SPH  windfield  and  Schloemer  pressure 
field  were  used  to  drive  the  hydrodynamics.  No  inflows  or  rainfall  were  con¬ 
sidered.  Simulated  water  level  histories  with  respect  to  an  initial  lake 
level  of  1-t.O  ft  msl  (plot  zero)  are  compared  with  observed  level  histories  in 
Plates  3-101  through  B-109.  Simulated  waterlevels  exceed  measured  water 
levels  in  the  south-western  section  of  the  lake  and  are  in  good  agreement  with 
observed  levels  in  the  other  lake  sections.  At  HGS-6  (Okeechobee)  the  simu¬ 
lated  peak  is  approximately  L.5  feet  less  than  the  observed  maximum.  Rainfall 
is  estimated  to  account  for  approximately  .3  foot  of  this  difference.  The 
timing  of  the  simulated  peaks  and  troughs  in  water  level  histories  is  in 
igreement  with  observed  data.  The  general  storm  track  and  windfield  structure 
,  r  -du.  o:  ’tic  <PH  simulation  is  thereby  substantiated.  Simulated  verticallv 
••'grtt-  a  urr.  .t  ire  shown  in  Plates  B  —  1  10  through  B-llb.  Global  current 
■  .  t  .re  :  present  oil  in  Plates  B— 117  through  B-1JA.  Maximum  current 

i  g:\  1 1  :  -»  >n  the  order  of  «  .  3  ips.  Simulated  significant  wave  heights  and 

,  •  r  :  -ire  mpareJ  with  measured  data  in  Table  A-J7.  Maximum  simulated  sig- 

.:  :  tit  wave  he  ignis  are  approximate  lv  <s  feet  with  periods  of  approximately 
•  Minimum  -.urge  levels  are  presented  i:.  Table  A-d8.  Boxed  cells 

:  :  r  ,,•••.•  ;  ;  .o,.  ;  ireas .  Areal  extent  of  simulated  flood  areas  corresponds  to 

:  i  •.  W-i'  "  ■  M  ]  .  Minimum  stilled  water  depths  are  shown  in  Table  A-J9 . 

:  :•  ;>r»  .  or  ireas  initially  wet,  which  dried  sometime  during 

.."..1  it  ;  .  A  minimum  depth  oi  1  foot  was  used  to  define  these  areas, 

i re  1 1  •  x'.-.-nt  I  exposed  areas  is  in  general  agreement  with  GW- 1  h 7  [31] 


■r.p  i  r  l  son  ■ !  i’31.  i:id  SPH 
over  •  :vn  imi  simulation  results 

J , j ,  SPH  simulated  water  levels  exceed  PBL  levels  at  all  gage  stations. 
'PH  r  *  *  s  ;  1 1  s  correspond  closer  to  measured  results  in  the  south  and  north- 
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eastern  sections  of  the  lake.  PBL  results  correspond  closer  to  measured 
results  in  the  southwestern  section  of  the  lake.  Wave  heights  and  periods  are 
larger  for  the  SPH  simulation  than  these  produced  in  the  PBL  simulation.  As 
mav  be  noted  by  comparing  Table  A-25  and  A-26  with  Tables  A-23  and  A-29, 
areas  of  inundation  and  exposure  are  larger  in  areal  extent  in  the  SPH  simula¬ 
tion.  All  of  these  results  are  a  direct  consequence  of  the  near  eye  structure 
of  the  two  windfields.  As  previously  noted,  the  near  eye  structure  is  quite 
different  and  SPH  wind  velocities  exceed  PBL  velocities  by  as  much  as  20  knots 
in  magnitude.  Directional  characteristics  are  generally  the  same,  however  may 
differ  by  up  to  25-30°. 

216.  It  is  a  difficult  task  in  general  to  simulate  the  overwater 
windfield  over  Lake  Okeechobee  for  the  August  1949  hurricane  due  to  eye 
pressure  adjustments  to  land  and  subsequent  readjustments  to  water  and  then 
readjustments  to  land.  It  appears  based  upon  these  two  simulations  that  the 
major  unknown  and  source  of  discrepancy  between  measured  and  simulated  water 
level  fluctuations  is  the  windfield  itself. 

247.  In  order  to  observe  the  magnitude  of  this  discrepancy  in  general, 
two  additional  storms  are  considered  in  turn  in  subsequent  parts. 


PART  VII:  OCTOBER  1950  HURRICANE 


C^-iier a  1  st..'  r m  oh a_r  ic  ter i st  ics 

2-8  .  lae  Oct 'her  17-18,  1950  (.Ring)  Hurricane  was  considered  to  further 
verifv  the  r c recast ing/h indcast ing  capability  of  the  complete  modeling  system. 
Hie  SPH  hurricane  model  formulation  was  used  to  produce  10-meter  10-minute 
average  sustained  overwater  windspeeds  to  drive  the  hydrodynamic  model.  The 
pvervater  windfieLds  were  reduced  over  Lake  Okeechobee  dvnamicallv  based  upon 
land/water  boundaries  developed  in  the  hydrodynamic  simulation. 

id 9.  The  October  17-18,  1950  (King)  Hurricane  possessed  a  maximum  cen¬ 
tral  pressure  depression  of  1.79  in  Hg  while  off  the  Florida  coast.  The 
radius  of  the  storm  was  onlv  6  nm  and  the  forward  speed  ranged  from  8.9  to 
11.7  kts.  The  storm  track  is  presented  in  Figure  VII-l.  As  may  be  observed, 
the  storm  passed  directly  over  Lake  Okeechobee  entering  east  of  Clewiston  and 
exiting  west  of  Okeechobee.  The  position  of  the  storm  is  given  in  Table  V 1 1  —  1 
in  hours  after  the  initial  position  cor respond ing  to  2300  HRS  Greenwich  Me  in 
Time  on  October  17,  1950. 

oPH  hurricane 

sub-modeL  simulation  resu Its 

250.  Twenty-four  hours  of  windfield  information  is  computed  esnpU  •-  ing 
three  snapshots.  Snapshot  data  employed  in  the  SPH  hurricane  simulation  are 
presented  in  Fable  VI 1—2 .  Snapshot  windfield  results  are  presented  in 
Plates  3-125  through  B  —  1 2  7  in  graphical  form.  The  entire  spatial  windfield  is 
presented  md  entered  with  respect  to  the  eye  in  these  figures.  Observe  the 
inf  1  w  angle  is  . enstant  at  each  radial  distance.  The  isovel  patterns  ire 
-'••.own  in  PI  ites  3-128  through  B-130  for  each  corresponding  snapshot.  The 
wind  era  t  ire  is  known  onlv  on  the  first  two  innermost  grids  of  the  live 
level  nest.  single  digit  numbers  represent  isovel  contours  at  a  Motor  of 

l.  Recalling  the  land/water  boundary  and  the  windfield  over  the  lake 
inter  i  ■  -  :  ->  i  ivnamic  feedback  svstom,  consider  results  for  the  windfiek! 

-o  1  -at  ;  md  pressure  field  at  gage  locations  shown  in  Figure  VI -2  as  pre- 
•  M  i’l  lies  3-1.31  through  B-199.  The  spatial  structure  of  the  winds 

r  Lake  I'kee  riobee  is  presented  in  Plates  B- 145  through  B-130.  Note  the 
sinus  in,  reuse  overwater  as  thev  feel  the  influence  of  the  lake.  Computed  and 
•■served  winds  ire  compared  in  Table  A- 30  for  snapshot  times  in  Table  VI I -2 . 
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October  19  50  (Day  1"  Hour  -3  CMT )  Position  I'ata 


T  itr.e 

0. 


Latitude  (oN) 
24.56 
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6. 

1 1 . 


12. 


25.3 
26.6" 
26.33 
2  "  .  0  ' 
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'omputed  wind  characteristics  are  in  general  agreement  v :  ?  •  ~e.as.re  a  ,e<  it 
cost  locations. 

-'PH  hydrodynamic 

sub-model  simulation  results 

2  5  2  .  A  period  of  23.5  hours  starting  at  1800  HRS  EST  on  1'  ctnher  l1' 
was  simulated  using  a  time  step  of  150  seconds.  Convective  acceleration  terms 


were  considered  in  the  middle  of  the  lake.  The  previously  presented  SPH  wind- 
field  and  the  Schloemer  pressure  field  were  used  to  provide  meteorologies  1 


•  .  v.  <•  .  n:  igur.it  ion  was 

:  •  :  .  1  iti-d  water  level 

,  ■  :  :  -.si  '  plat  area)  are 
.  -  i'.-i'l  through  3-163. 
wit:,  observed  levels 
■  •  r-  u  elewiston.  Simulated 
.  i'l  i t e s  3-lo4  through 
•. . -.  i : 1  •.  • .  1  tps  .  The  general 

.  :•  ■  .:  in  Plates  3-171  through 

t  •  the  corresponding  over  lake 
•  \r  .gh  3-130.  respectively.  The 
:  en t  structure.  Simulated 
pared  with  measured  data  in 
-ones  the  single  shoreline  wave 
are  less  than  1.0  ft  with  periods 
oding  occurred  during  the  October 
Islands  remained  dry. 

*0 
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PART  VIII:  SEPTEMBER  1979  HURRICANE 


General  storm  characteristics 

253.  The  September  3-4,  1979  (Hurricane  David)  was  considered  in  order 
to  further  test  the  hindcasting  capability  of  the  total  modeling  system  under 
the  present  levee  configuration.  The  SPH  hurricane  model  formulation  was  used 
again  to  produce  10-meter  10-minute  average  sustained  overwater  windspeeds. 

254.  The  September  3-4,  1979  (Hurricane  David)  possessed  a  maximum  cen¬ 
tral  depression  of  1.18  in  Hg  while  paralleling  the  east  coast  of  Florida. 

The  storm  radius  was  12  nm  with  a  forward  speed  of  10  knots.  The  storm  track 
is  presented  in  Figure  VIII-1.  The  storm  paralleled  the  east  coast  of  Florida 
eventually  making  landfall  near  Savannah,  Georgia.  The  position  of  the  storm 
is  given  in  Table  VIII-1  in  hours  after  the  initial  position  corresponding  to 
0400  HRS  Greenwich  Mean  Time  on  September  3,  1979. 

SPH  hurricane 

sub-model  simulation  results 

Thirty-six  hours  of  windfield  information  is  computed  using  one  snapshot, 
whose  characteristics  are  presented  in  Table  VIII-2.  Snapshot  windfield 
results  are  shown  in  Plate  B-  77  in  graphical  form,  with  constant  inflow 
angle  centered  about  the  eye.  The  corresponding  isovel  pattern  is  shown  in 
Plate  B-78  in  previous  format.  Single  digit  numbers  represent  isovel  con¬ 
tours  at  a  factor  of  ten  in  knots. 

255.  Noting  the  land/water  boundary  and  the  windfield  interact  as  a 
dyanmic  feedback  system,  consider  results  of  the  windfield  solution  and  pres¬ 
sure  field  at  gage  locations  shown  in  Figure  VI-2  as  presented  in  Plates 
B-79  through  B-92.  The  spatial  structure  of  the  overlake  winds  is 
presented  in  Plates  B-93  through  B-99.  Note  the  winds  increase  overwater. 
Computed  and  observed  winds  are  compared  at  Okeechobee  in  Table  A-32  for 
windfield  structure  times.  Computed  winds  are  somewhat  less  than  observed 
winds.  In  reference  to  Table  VIII-2,  for  a  central  pressure  deficit  of  40  mb 
the  SPH  winds  are  a  maximum  of  67  knots,  which  is  less  than  the  85  knot 
maximum  winds  shown  in  Figure  VIII-I  in  the  vicinity  of  the  lake. 


SPH  hydrodynamic 

sub-model  simulation  results 

256.  A  period  of  30  hours  starting  at  0400  HRS  GMT  on  3  September  1979 
was  simulated  using  a  time  step  of  300  seconds.  Convective  acceleration  terms 
were  considered  in  the  middle  of  the  lake.  The  previously  presented  SPH  wind- 
field  and  Schloemer  pressure  profile  were  used  to  provide  meteorological 
forcing.  No  inflows  or  rainfall  were  considered.  The  levee  configuration 
including  tree  islands  shown  in  Figure  V-3  was  considered.  Simulated  water 
levels  with  respect  to  an  initial  lake  level  of  14.48  msl  (plot  area1)  are  com¬ 
pared  with  observed  water  level  histories  in  Plates  B-200  through  B-203. 
Simulated  water  levels  are  in  reasonable  agreement  with  observed  levels  in 
overall  shape.  The  simulated  level  at  Port  Mayaca,  however,  does  not  reach 
the  level  observed.  This  is  due  to  lower  wind  loading  than  observed  over  the 
lake.  As  previously  noted,  for  a  40  mb  pressure  depression  one  would  not 
expect  an  85  knot  maximum  sustained  wind  magnitude.  It  would  appear  an 
adjustment  to  the  reported  control  pressure  depression  would  be  warranted. 
Simulated  vertically  integrated  current  histories  are  shown  in  Plates  B-204 
through  B-210.  Maximum  current  magnitudes  are  on  the  order  of  1.0  fps  thus 
facilitating  a  300  second  time  step.  The  general  current  structure  over  the 
entire  lake  is  presented  in  Plates  B— 2 1 1  through  B-217.  Simulated  maximum 
significant  wave  height  at  Clewiston  is  approximately  1.7  feet.  No  measured 
wave  data  are  available.  No  significant  flooding  occurred  during  Hurricane 
David,  although  Kreamer  an  Torrv  Islands  were  slightly  inundated  and  sections 
of  the  southern  group  to  Tree  Islands  were  overtopped  and  became  submerged. 

No  overtopping  of  the  major  levee  system  occurred. 
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PART  IX:  DESIGN  HURRICANES 


General  outline 

257.  The  development  of  both  standard  project  and  probable  r.aximum  hur¬ 
ricanes  is  presented  for  Lake  Okeechobee  in  terms  of  both  track  and  hurricane 
intensity  parameters.  A  single  standard  project  and  probable  maximum  hurri¬ 
cane  are  simulated  using  the  complete  hurricane  modeling  package;  i.e.,  the 
hurricane  sub-model  is  exercised  to  provide  input  for  the  hydrodynamic  model 
simulation.  In  the  final  sections  of  this  part,  the  effect  of  levee  breeching 
is  demonstrated  for  a  hypothetical  failure  near  Okeechobee  during  the  probable 
maximum  hurricane. 

Development  of 

design  hurricane  parameters 

258.  Schwerdt  et  al.  [ 21  have  developed  meteorological  criteria  for  the 

most  severe  hurricane  reasonably  characteristic  of  a  region,  Standard  Project 

Hurricane  (SPH) ,  and  for  the  hurricane  that  will  produce  the  highest  sustained 

wind  that  can  probably  occur  at  a  specified  coastal  location,  Probable  Maximum 

Hurricane  (PMH) .  A  single  limiting  value  for  the  meteorological  parameters  of 

peripheral  pressure  (P  )  and  central  pressure  (P  ) ,  while  upper  and  lower 

w  o 

limits  for  the  radius  to  maximum  winds  (R) ,  forward  speed  (T) ,  and  track 
direction  (0)  were  determined  at  100  nautical  mile  increments  as  shown  in  Fig¬ 
ure  IX— 1  along  the  Gulf  and  eastern  coasts  of  the  United  States. 

259.  The  central  pressures  adopted  for  the  SPH  and  PMH  are  shown  in 
Figures  IX-2  and  IX-3.  A  peripheral  pressure  (Py)  of  29.77  in.  (1008  mb)  was 
adopted  for  the  SPH,  while  a  peripheral  pressure  of  30.12  in.  (1020  mb)  was 
adopted  for  the  PMH.  The  adopted  upper  and  lower  limits  of  the  radius  to 
maximum  winds  are  shown  in  Figures  IX-4  and  IX-5  for  the  SPH  and  PMH,  respec¬ 
tively.  Adopted  SPH  and  PMH  upper  and  lower  limits  of  forward  speed  are  shown 
in  Figures  IX-6  and  IX-7,  respectively.  In  interpreting  the  maximum  allowable 
range  of  track  direction  for  the  SPH  and  PMH  shown  in  Figures  IX-8  and  IX-9, 
respectively,  the  relations  between  direction  and  forward  speed  shown  in 
Table  IX-l  must  be  used.  The  track  direction  represents  the  angle  from  which 
the  storm  is  approaching  the  coast  measured  clockwise  from  North;  i.e.,  a 
meteorological  definition  is  employed. 

260.  In  developing  SPH  and  PMH  storms  for  Lake  Okeechobee,  the  condi¬ 
tions  at  the  appropriate  coastal  sections  are  extrapolated  to  the  vicinity  of 
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For  the  PMH 


Speed  cateeor% 


For  the  SPH 


Speed  category 


Forward  speeds  (T) 

6  kt  <  T  _<  10  kt 
(11  km/hr  <_  T  £  19  km/hr) 

10  kt  <  T  <_  36  kt 
(19  km/hr  <  T  <_  67  km/hr) 

T  >  36  kt 
(T  >  67  km/hr) 


Forward  speeds  (T) 


4  kt  <  T  <  10  kt 
(7  km/hr  <_  T  _<  19  km/hr) 

10  kt  <  T  <  36  kt 
(19  km/hr  <  T  <_  67  km/hr) 

T  >  36  kt 
(T  >  67  km/hr) 


the  lake.  Design  geometry  and  coastal  design  storm  characteristics  are  pre¬ 
sented  prior  to  developing  these  general  extrapolation  procedures. 

Design  geometry 

261.  In  order  to  specify  the  orientation  of  design  hurricanes  relative 
to  the  lake,  it  is  necessary  to  specify  certain  geometrical  relations.  The 
latitude  and  longitude  ordered  pair  (A,  L)  of  the  center  of  the  lake  is  speci 
fled  in  conjunction  with  a  critical  design  radius,  r  .  The  critical  design 
radius  represents  the  approximate  distance  in  nautical  miles  of  the  average 
hurricane  radius  to  maximum  winds  known  to  have  influenced  the  lake.  From 
previous  work  r  20  nm  ,  A  »  26.92°N  ,  L  -  80.83°W  . 

262.  In  general,  it  would  appear  necessary  to  consider  storms  passing 
through  each  point  on  the  circle  shown  in  Figure  IX-10  from  all  compass  direc 
tions.  By  consulting  the  limits  of  track  direction  for  coastal  segments  for 
the  SPH  and  PMH  given  by  Schwerdt  et  al.  [2],  it  is  possible  to  eliminate 
certain  angle  bands. 
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Figure  IX-12.  Coastal  Segment  Determination 
of  Design  Hurricanes 


263.  If  one  assumes  that  the  most  critical  storm  passes  through  the 
given  point  tangent  to  the  circle  at  that  point,  it  is  possible  to  develop  the 
design  condition  shown  in  Figure  IX-11.  For  each  point  on  the  circle  with 
critical  design  radius  there  is  one  critical  storm  direction.  In  addition, 
all  points  on  the  circle  need  not  be  considered  based  upon  the  limits  of  track 
direction.  The  limits  on  track  direction  define  the  locus  of  points  on  the 
circle,  design  circle  sector,  between  140-340°  clockwise  from  North.  Under 
this  convention,  the  lake  center  is  located  always  one  critical  design  radius, 
r  ,  to  the  right  of  the  storm  track.  This  orientation  is  assumed  critical. 


Coastal  design  storm  characteristics 


264.  The  coastal  points  shown  in  Table  IX-2  encompass  Lake  Okeechobee 
as  shown  in  Figure  IX-1.  The  latitude  and  longitude  coordinates  of  these 
points  are  as  given  in  Table  IX-2.  The  SPH  and  PMH  desgin  characteristics  are 
shown  in  Tables  IX-3  and  IX-4,  respectively  for  the  coastal  segments  through 
which  storms  pass  in  influencing  Lake  Okeechobee. 


General  extrapolation  procedures 


265.  Given  the  latitude  and  longitude  (X,  L)  of  the  center  of  the  Lake, 
the  critical  design  radius,  r  ,  and  limits  of  the  design  circle  sector  angles 
and  ®nax  measured  as  compass  bearings,  one  considers  arbitrary  points 


Coastal 


Table  IX-2 


segment  Points 


Location 


Coastal  Distance  (nm) 


Latitude  (°N’) 


Long 


Table  IX-3 

SPH  Design  Characteristics 


Segment  No. 


Coastal  Distance  Ranee 


1200-1300 

1300-1400 

1400-1500 

1500-1600 


27.42 

27.19 

27.12 

27.25 

17.08 


30.5  4. 

29. 

28.5 

30.  4. 


29.77  in. 


Table  IX-4 

PMH  Design  Characteristics 


W.*.'  i*!1 


where 


:'p  =  Latitude  of  point  P 
3  =  Latitude  of  center  of  the  lake 
Lp  =  Longitude  of  point  P 
L  =  Longitude  of  center  of  the  lake 

a  =  Compass  bearing  of  point  P  (clockwise  from  North) 
r  =  Design  radius  in  nautical  miles 
(Note:  1'  of  arc  on  the  earth  corresponds  to  1  nautical  mile) 

266.  One  next  constructs  the  line  passing  through  point  1  tangent  to 
the  circle  of  that  point  as  shown  in  Figure  IX-12.  The  general  formula  for 
the  equation  of  this  line  in  latitude  and  longitude  space  is  given  as  follows: 


mD 


x  +  b„ 


where 

Bp  >  -  tana 
bD  =  Lp  mDAp 

One  next  considers  each  coastal  segment  i  »  i  ■  l,. ..4  ,  as  shown  in  Fig¬ 
ure  IX-12.  Each  segment  is  represented  as  a  straight  line  of  the  form  given 
below. 


y  =  nijX  +  b 


where 


mi  =  (yi+L  "  yi)/(xi+l  “  xi} 


b  =  y  -  mb 
i  i  i  x 

The  line  representing  the  storm  is  next  checked  against  each  coastal  segment 
i  until  an  intersection  with  segment  i  is  found.  The  procedures  are  out¬ 
lined  as  follows: 

1.  If  m^  =  m^  the  lines  are  parallel  and  do  not  intersect.  Skip  to 
segment  i  +  1  . 

2.  For  m,  j4  ra  x,  -  b,  -  b  /b,  -  m  compute: 

i  o  1  i  o  i  o  r 

3.  If  x  £  x  ^  x.  ,  ,  the  lines  intersect.  If  x  <  x.  or  x  >  x.  . 

1  I'M  1 

skip  to  segment  i  +  1  . 


ss 


4.  If  the  lines  intersect: 

a.  Employ  (print)  the  SPH  and  PMH  design  characteristics  for 
coastal  segment  i  . 

b.  Compute  y^  "  m^ 

c.  Compute  the  distance  between  the  point  of  intersection  I  and 
point  P  as  follows: 


(X  -  Xr)2  +  [ (L  -  Yr)  cos  Xr ] 2 


d.  Next  calculate  the  time  since  landfall  that  it  takes  the  storm 
center  to  reach  point  P  : 


e. 

f . 

tr¬ 

iable  IX- 5 

Adjustment  Factors  For  Reducing  Hurricane  Wind  Speeds 
When  Center  Is  Overland  (Category  B) 


Hrs  After  Reduction 

Center  Is  Overland  Factor 

7.5  .90 

15.0  .80 

22.7  .70 


X  =  — 

P  Ti 


(Note  T  is  in  hrs,  since  d 


is  in  (nm)  and  is  in  kts) 


Enter  Table  IX-5  and  determine  and  print  the  reduction  factor 
and  T 

P 

Proceed  to  the  next  point  P  =  8  +  A0 

P 


Program  DESIGN 

267.  The  above  algorithm  is  employed  in  Program  DESIGN  to  develop 
design  storm  characteristics  for  the  lake.  A  listing  of  program  along  with 
attendant  job  control  language  and  input  data  is  presented  in  Appendix  E. 
,»\V  Final  design  storm  characteristics  are  presented  for  angle  a  -  140°  in 


Table  A-33. 


Development  of  a 
Standard  Project  Hurricane 

168.  Based  upon  characteristics  associated  with  coastal  segment  1300- 
1400,  a  hurricane  was  specified  as  shown  in  Table  A-34 .  The  joint  probability 
method  storm  track  geometry  described  in  Part  III  of  the  second  report  of  this 
series  was  employed  in  the  hurricane  sub-model.  A  meteorological  angle  of  90° 
corresponding  to  a  compass  bearing  of  180°  was  used  to  specify  the  track  of  a 
storm  passing  through  the  center  of  the  lake.  Nine  hours  of  winds  were  devel¬ 
oped  using  three  snapshots  in  previous  format  as  shown  in  Plates  B-192 
through  B-197  by  the  parametric  standard  project  hurricane  formulation.  The 
initial  storm  position  is  along  the  track  at  a  distance  from  lake  center  equal 
to  the  maximum  of  three  times  the  radius  of  maximum  winds  and  120  nautical 
miles.  The  final  storm  position  is  along  the  track  at  a  distance  from  lake 
center  of  three  times  the  radius  of  maximum  winds  in  nautical  miles.  Snapshot 
windfield  vector  and  corresponding  isovel  contours  are  presented  in  Plates 
B-218  through  B-221  using  previously  described  formats.  Maximum  winds  exceed 
100  knots.  Observe  ^n  Plates  B-222  and  B-223  that  the  central  pressure 
deficit  has  decreased  from  2.58  in  Hg  to  2.01  in  Hg  due  to  the  hurricane 
filling  relationships  in  Table  IX-5. 


Standard  project  hurricane 
initial  hydrodynamic  simulation  results 

269.  In  order  to  study  the  hydrodynamic  response  of  the  lake  to  the 
standard  project  hurricane  previously  developed,  a  numerical  simulation  was 
performed  with  the  lake  at  an  initial  level  of  15.5  msl.  Four  hundred  thirty 
time  steps  of  75  seconds  duration  were  used  to  simulate  8.95  hours.  No 
inflows  or  rainfall  were  considered.  The  present  levee  and  tree  island  con¬ 
figuration  shown  in  Figure  V-3  was  considered.  In  this  initial  simulation, 
windfield  and  pressure  information  were  updated  every  24  time  steps  corre¬ 
sponding  to  0.5  hour.  Simulated  water  level  histories  are  presented  with 
respect  to  15.5  msl  plot  zero  in  Plates  B-224  through  B-230  for  HGS-1  through 
HGS-6  and  Port  Mayaca,  respectively.  Peak  water  levels  with  respect  to  mean 
sea  level  and  their  associated  time  of  occurrence  are  presented  in  Tables  A-35 
and  '-^b,  respectively.  Storm  loadings  are  presented  in  Plates  B— 2 3 1  through 
B-244  .  Observe  the  stair  step  type  description  at  half  hour  intervals.  The 
current  response  due  to  the  spatial  distribution  of  overlake  winds  is  por¬ 
trayed  in  Plates  B-245  through  B-260  at  times  of  l,  2,  3,  4,  5,  6,  7,  and 


3.96  hours,  respectively.  During  the  first  four  hours  of  the  storm,  winds  are 
directed  toward  the  western  side  of  the  lake.  The  water  surface  is  elevated 
at  the  west  and  lowered  on  the  east  side  of  the  lake.  Water  tends  to  pile  up 
against  the  southern  tree  islands  from  the  shore  side  and  some  overtopping 
occurs  from  the  shore  side.  For  the  northern  band  of  tree  islands,  over¬ 
topping  occurs  from  the  lakeside.  This  pattern  becomes  completely  reversed  as 
the  :  torm  passes  through  the  center  of  the  lake  and  the  winds  shift  toward  the 

east.  The  role  of  the  tree  islands  in  the  overall  response  of  the  lake  is 
very  complex.  As  may  be  noted  in  Plate  B-230  at  Port  Mayaca  the  resurgence 
at  time  7.8  and  8.4  hours  may  be  due  to  long-wave  reflection  from  the 
southeastern  entrance  of  the  tree  islands.  Maximum  currents  are  cn  the  order 
of  7.5  fps.  The  maximum  significant  wave  height  at  Clewiston  is  5.2  ft. 


Standard  project  hurricane 
additional  hydrodynamic  simulations 

270.  Two  additional  hydrodynamic  simulations  were  performed  using  all 
the  conditions  of  the  initial  simulation.  In  simulation  one,  the  wind  and 
pressure  fields  were  interpolated  at  a  variable  number  of  time  steps  based 
upon  the  location  and  size  of  the  storm.  These  procedures  are  more  fully 
developed  in  Part  X.  Water  surface  histories  with  respect  to  15.5  msl  plot 
zero  and  the  corresponding  wind/pressure  histories  are  given  in  Plates  B -  2 6 1 
through  B-266  at  Clewiston,  Canal  Point,  and  Okeechobee,  respectively.  In 
simulation  two,  the  wind  and  pressure  fields  are  interpolated  every  time  step. 
Water  surface  histories  and  corresponding  wind/pressure  histories  are  given  in 
Plates  B-267  through  B-  272.  Water  surfaces  histories  for  each  additional 
simulation  are  equivalent  to  the  initial  simulation  results.  Computational 
time  requirements  are  vastly  different  as  shown  in  Table  A-37. 

271.  It  would  appear  the  approach  used  in  the  initial  simulation  or  in 
additional  simulation  one  is  considerably  cheaper. 


Development  of  a 
probable  maximum  hurricane 

272.  Based  upon  characteristics  associated  with  coastal  segment  1  300- 
1400  nm,  the  hurricane  developed  in  Table  A-38  is  developed.  A  meteorological 
angle  of  90°  corresponding  to  a  compass  bearing  of  180°  was  used  to 
specify  the  track  of  a  storm  passing  a  distance  of  one  radius  to  maximum  winds 
to  the  left  of  the  lake  center.  Nine  hours  of  winds  were  developed  using 
three  snapshots  in  previous  format  as  shown  in  Plates  B  —  2 7 3  through  B-  278  by 


the  parametric  standard  project  hurricane  formulation.  The  initial  storm 
position  is  along  the  track  at  a  distance  of  120  nautical  miles  from  the 
meteorological  line.  The  final  storm  position  is  along  the  track  at  a  dis¬ 
tance  from  the  meteorological  line  of  three  times  the  radius  of  maximum  winds. 
Note  maximum  sustained  windspeeds  are  greater  than  130  knots. 

Probable  maximum  hurricane 

initial  hydrodynamic  simulation  results 

273.  In  order  to  investigate  the  hydrodynamic  response  of  the  lake  to  a 
hurricane  of  probable  maximum  intensity,  a  numerical  simulation  was  performed 
with  the  lake  at  an  initial  level  of  24.5  msl  corresponding  to  the  standard 
project  flood  elevation.  A  45  second  time  step  was  employed  in  simulating 
3.5  hours.  No  inflows  or  rainfall  were  considered.  The  present  levee  and 
tree  island  configuration  shown  in  Figure  V-3  was  considered.  In  this  initial 
simulation,  windfield  and  pressure  information  were  updated  every  40  time 
steps  corresponding  to  0.5  hour.  Simulated  water  level  histories  with  respect 
to  24.5  msl  plot  zero  are  presented  in  Plates  279-285  for  HGS-1  through  HGS-6 
and  Port  Mayaca,  respectively.  Peak  water  levels  with  respect  to  mean  sea 
level  and  their  associated  time  of  occurrence  are  presented  in  Tables  IX— 1 J 
and  IX-13,  respectively.  Storm  loadings  are  presented  in  Plates  B- 286 
through  S-299  in  a  stair  step  manner  at  half  hour  intervals.  The  current 
response  due  to  the  spatial  distribution  of  overlake  winds  is  portrayed  in 
Plates  B-300  through  B-3I7  times  of  1,  2,  3,  4,  5,  6,  7,  8,  and  8.5  hours, 
respectively.  The  response  from  hour  7  on  is  such  that  lake  currents  are 
opposite  in  direction  to  wind  loading  over  the  eastern  half  of  the  lake. 

Phis  is  apparently  due  to  the  fact  that  the  wind  loading  is  not  sufficient  to 
overcome  the  water  surface  elevation  gradient  developed  during  the  previous 
hours  of  the  storm.  At  an  initial  lake  level  of  24.5  msl  all  tree  islands  are 
initially  submerged.  As  water  tends  to  drawdown  faster  on  the  lake  side  of 
the  southern  tree  islands,  some  sections  are  overtopped  from  the  shore  side 
during  the  early  hours  of  the  storm.  Eventually  after  4-5  hours  all 
southern  tree  islands  become  exposed.  Maximum  lake  currents  are  on  the  order 
to  12.3  ips.  The  maximum  significant  wave  height  at  Okeechobee  is  9.97  feet. 
At  hour  seven,  the  levee  at  Okeechobee  is  initially  overtopped.  The  still 
water  level  eLevation  exceeds  the  levee  height  from  hour  7  through  hour  8. 
During  this  period,  the  overtopping  contribution  due  to  waves  is  in  the  range 
of  -*-6  1  cubic  feet  per  second  per  foot  of  levee. 


Probable  maximum  hurricane  levee 
breeching  hydrodynamic  simulation 

27-t.  7he  section  of  levee  at  Okeechobee  was  specified  to  be  in  i 
breeched  condition  at  hour  7  based  upon  the  initial  probable  maximum  hurri  anc 
results.  This  initial  simulation  was  rerun  with  a  levee  breech  at  hour  7  ; 

duration  of  six  hours.  The  breech  was  specified  near  Okeechobee  as  shown  in 
Figure  V-3  as  shown  in  Table  A-4 1 .  The  hydrodynamic  response  of  the  lake 
remained  essentially  unchanged  as  may  be  noted  by  comparing  Table  A— *2  and 
A--*3  with  Table  A-39  and  A-40  respectively.  The  peak  water  level  in  cell 
1 33 ,  3)  at  Okeechobee  is  reduced  by  .1  feet  only.  Representative  breech 
conditions  are  shown  in  Table  A— id  at  hour  6.8875  (before  the  breech),  at 
hour  7.13  (immediately  after  the  breech),  and  at  hour  8.023.  Note  water 
levels  are  sufficient  to  overtop  both  barriers  prior  to  the  initiation  of  the 
nreech.  Immediately  after  the  breech,  flow  over  section  (38,  3)  reaches  over 
300,000  cfs,  while  at  section  (39,  3)  the  flow  is  approximately  60,000  cfs. 

By  hour  8.023,  the  water  level  at  the  northern  end  of  the  lake  has  lowered 
such  that  only  section  (38,  3)  is  overtopped  at  a  rate  of  approximately  23,000 


: :  s .  Beyond  hour  8.162u,  no  overtopping  occurs.  Levee  elevations,  under  the 
urrent  model  formulation,  continue  to  be  lowered  at  a  uniform  rate  per  time 
step,  determined  by  subtracting  the  initial  elevations  in  Table  IX-4  from 
the  final  elevations  and  dividing  the  result  by  the  duration  of  breech  then 
-ult  inlying  the  time  step  length. 


I  >  i t  nr oba_b  i_l  i ty  method  considerations 

273.  As  developed  in  Report  2  in  this  series,  approximately  600  com- 
:  in.ed  hurricane  sub-model  hvdrodvnamic  sub-model  simulation  pairs  must  be 
executed  in  i  joint  probability  method  design  for  present  levee  and  hurricane 
cite  structures.  In  the  previous  sections  of  Part  IX  a  variable  wind/pressun 
:ie!  :  interpolation  procedure  was  developed  producing  essentially  the  same 
results  ,3  the  single  time  step  procedures  used  in  Parts  V I  —V III  during  the 


.  :  r  it  i. 


process  but  at  substant ia i 1 v  less  cost.  As  an  additional  ccst- 


'.tving  or  .  dure,  it  is  recommended  that  a  variable  time  step  approach  be 
invest  ig ued .  This  approach  could  be  implemented  in  the  same  manner  as  the 
v  iriable  interpolation  procedure;  e.g.,  only  when  the  hurricane  is  close  to 
the  lake  is  the  time  step  reduced  from  a  user  specified  baseline  value 
i  yt>>  seconds)  .  The  entire  concept  of  storm  track  is  in  the  concept  M  the 


The  entire 


:oint  prababil  ity  net  hod  further  investigated  in  this  part. 

-7o .  in  Part  IX  paragraph  269,  a  standard  proiu  t  iiurri.  ui.-  *  is  >imul  jted 
employing  t  i  •  e  characteristics  shown  in  Table  A- i-  .  ::  i  ne  hours  correspond inc 

to  ISO  nautical  miles  of  track  starting  .it  i  position  Id?  nautical  miles  due 
south  from  lake  center  were  considered.  In  examining  Table  A- 3b ,  for  lake 
areas  near  HCS-3,  there  was  some  question  as  to  whether  or  not  peak  water 
levels  had  been  reached.  The  storm  track  algorithm  was  modified  is  shown  in. 
the  relation  below  to  allow  for  seicning  resurgence  effects. 

T  =  [  ma  k(120.,3R)+  20.  +  3  R  1  /  V  +  2.3  (9.1) 

where 

T  =  Total  storm  track  simulation  time  (hrs) 

R  *  Storm  radius  (nut) 

Vf  =  Storm  forward  speed  (kts) 

The  Lasc  term  representing  2.5  hours  represents  the  travel  time  of  i  long-wave 
from  one  side  of  the  lake  to  the  other.  A  generalized  hurricane  sub-model 
simulation  was  performed  using  the  above  refined  track  requirement  as  shown  in 
Table  A- 2  5 .  Windfield  information  is  developed  at  hours  Zi-ro  and  on..  It  is  the 
same  developed  previously  in  Plates  B- 218  and  B—  219  >nd  plates  B-  220  and 
3-221,  respectively. 

"7.  The  wind  information  at  hour  eleven  is  Vo  own  in  plates  B- 318  and 
3- 319,  respectively,  using  previous  formats.  Observe  the  CPI  has  been 
reduced  from  the  hour  nine  value  of  2.01  to  1.89  in  Hz  at  hour  eleven.  A 
corresponding  hydrodynamic  model  simulation  was  performed  for  i  period  of 
10.  a  hours.  Maximum  surge  elevation  (msl)  and  their  >.  or  responding  times  of 
occurrence  are  shown  in  Tables  A-26  and  A— 27 ,  respective!  v  .  in.  .mp.irinz 
these  results  with  those  shown  previously  in  Tables  A- 5  3  o'.d  .\-  ib ,  respec¬ 
tively,  very  little  differences  is  noted  even  in  the  vicir.it;.  a  HiCS-i.  How¬ 
ever,  if  one  observes  Plates  B-229  and  the  corresponding  r. -suits  shown  in 
Plate  B- 3J0 ,  it  is  necessary  to  simulate  over  nine  nours  in  order  to  be  -ure 
a  peak  level  has  been  reached  at  HGS-6 .  For  this  reason  end  t-  provide  i 
margin  of  safety  the  track  algorithm  given  in  Bquat  ion  1  1 . i  ,  ls  recommended. 


PART  X:  STUDY  RESULTS  AND  CONCLUSIONS 


~  ’  ^  ‘  i.ie.ii.  studv  results  are  presented  and  the  following  conclusions 
are  drawn: 

a.  A  generalized  modeling  system  has  been  developed  for  deter¬ 
mining  seiche  and  hurricane  induced  water  level  fluctuations 
in  Lake  Okeechobee.  The  system  consists  of  two  major  compo¬ 
nents:  l)  a  hurricane  sub-model  (Program  LOGHM)  and  a  attend¬ 

ant  graphics  software  (Program  HGRAPH) ,  and  2)  a  hydrodynamic 
submodel  description  of  a  long-wave  and  short-wave  phenomena 
(Program  LAKE)  complete  with  graphics  (Procedure  POST) . 

1).  The  hurricane  sub-model  consists  of  a  direct  parametric 

approach  as  well  as  a  planetary  boundary  layer  formulation. 
Both  approaches  are  completely  developed  in  this  report. 
Generalized  hurricane  track  interpolation  procedures  have  been 
developed  for  hindcasting/forecasting  as  well  as  joint  proba¬ 
bility  method  (straight  line)  storm  tracks. 

c.  The  hydrodynamic  sub-model  long  and  short  wave  theory  are  com¬ 
pletely  developed  in  this  report.  Wave  runup  and  overtopping 
as  well  as  levee  breeching  are  considered.  Within  the  long¬ 
wave  calculations  sub-grid  scale  barriers  and  a  moving 
boundary  are  treated. 

d.  A  global  grid  (64  x  48)  of  3072  cells  has  been  constructed  for 
hydrodynamic  sub-model  calculations.  SAJ  personnel  assisted 
in  developing  lake  topography  and  levee  profile  data  suitable 
for  the  global  grid. 

e.  The  August  12-13,  1968  seiche  and  the  August  20-21,  1964 
seiche  were  used  to  calibrate  and  verify  the  long-wave  bottom 
friction  mechanics  independent  from  wind  surface  loading. 
Simulated  levels  are  in  reasonable  agreement  with  measured 
water  surface  elevations  at  most  gaging  stations. 

Based  upon  simulations  results,  the  bottom  friction  mechanics 
including  canopy  effects  seem  reasonable.  In  the  seiche  simu¬ 
lations,  the  initial  lake  surface  elevation  is  accurately 
known  at  only  a  few  stations.  This  makes  the  specification  of 
the  initial  lake  profile  difficult.  Due  to  the  small  ampli¬ 
tudes  of  the  seiches  (approximately  .5  feet),  the  specifica¬ 
tion  of  the  initial  lake  surface  is  very  critical  for 
simulating  these  events. 

f.  The  August  1949  hurricane  was  considered  as  the  most  critical 
event  for  surface  wind  loading  calibration.  Both  the  PBL  and 
SPH  approaches  in  the  generalized  hurricane  sub-model  were 
exercised.  Although  peak  maximum  sustained  windspeeds  were  on 
the  order  of  3-8  knots  lower  in  the  PBL  approach,  the  inner 
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structure  (winds  within  a  vicinity  of  2  radii  of  maximum  winds 
from  the  storm  center)  was  quite  different.  Overlake  simu¬ 
lated  windspeeds  were  compared  with  measured  winds  at  several 
locations  around  the  lake.  Simulated  PBL  and  SPH  windspeeds 
were  in  excellent  agreement  with  measured  results  during  the 
initial  and  final  hours  of  th  storm  track.  During  the  period 
the  storm  center  was  near  or  over  the  lake,  PBL  windspeeds 
were  as  much  as  25  knots  lower  than  SPH  windspeeds  at  some 
gaging  stations.  SPH  windspeeds  were  in  reasonable  agreement 
with  observed  windspeeds  during  the  storm  period.  Simulated 
PBL  wind  directions  were  usually  within  30  degrees  of  simu¬ 
lated  SPH  wind  directions  for  the  gaging  stations  considered. 
Simulated  SPH  wind  directions  agreed  within  20-30  degrees  of 
observed  directions. 

The  hydrodynamic  sub-model  was  used  to  simulate  hurricane 
induced  water  levels  for  both  windfield  formulations.  Simu¬ 
lated  PBL  induced  water  levels  were  lower  than  SPH  induced 
water  levels  at  all  gaging  stations.  However,  PBL  generated 
levels  appeared  to  be  in  closer  agreement  with  observed  levels 
in  the  southern  portions  of  the  lake,  while  SPH  generated 
levels  were  in  closer  agreement  in  the  eastern  and  northern 
portions  of  the  lake.  The  maximum  water  level  excursion  at 
Okeechobee  was  approximately  14.5  feet.  The  PBL  generated 
maximum  water  level  excursion  at  Okeechobee  was  10.0  feet 
(31%  error),  while  the  SPH  generated  maximum  excursion  at 
Okeechobee  was  12.8  feet  (12%  error).  Simulated  significant 
wave  heights  were  slightly  larger  for  the  SPH  windfield  than 
for  the  PBL  windfield.  SPH  windfield  generated  significant 
wave  heights  were  in  reasonable  agreement  with  measured 
results.  Significant  wave  heights  and  periods  are  larger  for 
the  CW-167  curve  approach  than  under  the  SPM  methodology. 
Considering  the  CW-167  and  SPM  results  to  define  a  range, 
measured  significant  wave  heights  and  periods  fall  within  the 
range.  Regions  of  inundation  generated  by  the  SPH  windfield 
are  in  excellent  agreement  with  observed  regions.  For  the  PBL 
simulations,  regions  of  inundation  are  reduced  considerably  in 
areal  extent. 

Based  upon  the  results  produced  for  the  August  1949  hurricane, 
it  appears  that  the  wind  and  bottom  stress  formulations  are 
reasonable  in  the  long  wave  calculations.  In  order  to  further 
reduce  the  maximum  error  between  simulated  and  observed  water 
levels,  it  is  necessary  to  better  define  the  inner  structure 
of  the  hurricane  windfield  during  passage  over  the  lake.  The 
SPH  formulation  is  superior  to  the  PBL  formulation  based  upon 
hydrodynamic  simulation  results. 

The  October  1950  storm  although  not  of  the  intensity  of  the 
August  1949  hurricane  passed  directly  over  the  lake.  The  SPH 
approach  was  employed  to  generate  the  windfield.  Simulated 
far  field  windspeeds  and  directions  are  in  good  agreement  with 
measured  data.  Simulated  near  field  windspeeds  at  the  time  of 
lake  passage  are  within  20  knots  of  observed  data  in  all  cases 
and  at  most  stations  are  within  10  knots.  SPH  windfield 
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induced  peak  water  levels  produced  by  running  the  hydrodynamic 
sub-model  with  the  SPH  windfield  as  input  are  within 
I. -1.5  feet  of  peak  measured  levels  except  at  HGS-3.  Very 
little  measured  wave  data  are  available  at  shoreline  stations. 
Simulated  wave  data  are  in  agreement  with  the  limited  measured 
data . 


Based  upon  the  results  produced  for  the  October  1950  hurri¬ 
cane,  the  wind  and  bottom  stress  formulations  in  the  long  wave 
calculation  are  confirmed.  Discrepancies  between  observed  and 
simulated  water  level  histories  are  due  to  differences  between 
the  observed  and  computed  windfields.  Short  wave  simulated 
data  appear  reasonable. 


The  September  1979  hurricane  (1  avid)  was  considered  in  order 
to  check  the  model  under  the  piesent  levee  configuration. 
Hurricane  David  produced  overlake  winds  on  the  order  of 
60  knots;  i.e.,  winds  were  not  of  hurricane  intensity  over  the 
lake  unlike  the  August  1949  and  October  1950  storms.  The  SPH 
approach  was  used  to  generate  the  windfield  input  to  the 
hydrodynamic  sub-model.  Peak  simulated  overlake  winds  were 
40  knots,  approximately  20  knots  less  than  observed  values. 

In  the  vicinity  of  Lake  Okeechobee,  a  measured  central  pres¬ 
sure  depression  of  40  mb  is  reported  to  have  produced  maximum 
sustained  winds  of  85  knots.  Based  upon  normal  pressure  rela¬ 
tionships  for  calculating  maximum  overwater  winds,  an  85  knot 
peak  wind  cannot  be  generated.  The  SPH  windfield  maximum  was 
approximately  65  knots  with  overlake  maximum  windspeeds  at 
40  knots.  As  one  might  suspect,  SPH  generated  water  levels 
are  lower  than  observed  water  levels  for  gage  stations  nearest 
the  passage  of  the  storm;  i.e..  Port  Mayaca.  No  measured  wave 
data  are  available. 


It  is  difficult  to  draw  any  other  conclusion  based  upon  these 
results,  than  the  fact  that  the  hurricane  sub-model  input  data 
must  be  consistent  in  order  to  produce  reasonable  windfield 
results.  Peak  simulated  water  levels  are  directly  dependent 
on  the  hurricane  sub-model  windfield.  In  effect,  a  chain 
reaction  occurs  in  which  errors  in  one  step  of  the  process 
generate  errors  in  the  next  step. 


In  order  to  develop  evacuation  planning  alternatives,  a  design 
hurricane  methodology  was  developed  in  conjunction  with  the 
hurricane-hydrodynamic  model  pair.  Standard  Project  Hurri¬ 
canes  and  Probable  Maximum  Hurricanes  are  developed  for  Lake 
Okeechobee  by  formally  extrapolating  in  Program  DESIGN 
(Appendix  E)  design  hurricane  characteristics  developed  for 
Atlantic  and  Gulf  Coast  sections  by  Schwerdt  et  al.  [2].  A 
general  design  storm  radius  concept  is  developed  in  order  to 
specify  design  storms  passing  within  an  arbitrary  distance 
from  lake  center. 


A  Standard  Project  Hurricane  passing  directly  on  a  line  from 
South  to  North  through  the  center  of  the  lake  was  considered 
for  an  initial  level  of  15.5  msl.  The  SPH  approach  was  used 
to  generate  the  design  windfield.  Peak  sustained  winds 


exceeded  100  knots  on  the  lake.  A  hydrodynamic  simulation  was 
performed  using  the  SPH  generated  design  windfield.  Since  the 
initial  lake  level  corresponded  to  the  elevation  of  the  tree 
islands,  at  alternate  portions  of  the  simulation  sections  of 
the  tree  islands  were  overtopped,  exposed,  and  completely  sub¬ 
merged.  Hydrodynamic  sub-model  output  formats  were  extended 
to  include  maximum  and  time  of  maximum  water  levels  as  well  as 
corresponding  depth  fields.  A  variable  interpolation  interval 
scheme  was  developed  in  order  to  reduce  computational  cost  in 
the  hydrodynamic  sub-model.  Simulated  peak  water  levels 
exceeded  10.0  ft  above  the  12.5  msl  model  datum.  No  over¬ 
topping  of  the  levee  system  occurred. 

A  Probable  Maximum  Hurricane  passing  a  distance  of  one  radius 
to  maximum  winds  to  the  west  of  lake  center  on  a  direct  South 
to  North  track  was  considered  for  an  initial  lake  level  of 
24.5  msl  corresponding  to  the  Standard  Project  Flood.  Wind- 
speeds  in  excess  of  120  knots  are  developed  over  the  lake  by 
the  SPH  windfield  formulation.  This  storm  represents  the  worst 
possible  condition  for  Lake  Okeechobee.  Water  level  eleva¬ 
tions  of  over  19  feet  above  the  initial  lake  level  are  pro¬ 
duced  at  HGS-6  during  the  hydrodynamic  sub-model  simulation. 
Maximum  significant  wave  heights  on  the  order  of  10  feet  are 
produced  at  HGS-6.  The  levee  at  HGS-6  is  overtopped  for  a 
period  of  approximately  one  hour. 

Based  upon  the  simulation  of  this  event,  the  complete 
hurricane-hydrodynamic  model  pair  is  capable  of  simulating  all 
storm  events  associated  with  a  joint  probability  analysis. 

The  Probable  Maximum  Hurricane  condition  just  reported  was 
reconsidered  for  the  specification  of  a  levee  breech  at 
Okeechobee.  The  breech  was  specified  to  occur  at  hour  seven 
and  last  for  the  next  six  hours  as  shown  in  Table  1X-4 .  The 
hydrodynamic  response  of  the  lake  remained  essentially 
unchanged.  The  peak  water  level  at  Okeechobee  was  reduced  by 
.  1  feet . 

This  simulation  demonstrated  that  levee  breech  flows  may  be 
calculated  as  input  to  a  link-node  drainage  model  for  flood 
routing  and  emergency  action  planning. 

The  specifications  of  the  track  for  joint  probability  method 
application  was  refined  to  allow  the  travel  of  a  long-wave 
across  the  entire  lake  length  in  order  to  study  resurgence 
effects  on  peak  water  levels. 
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